Wide-field fluorescence molecular tomography with compressive sensing based preconditioning

: Wide-field optical tomography based on structured light illumination and detection strategies enables efficient tomographic imaging of large tissues at very fast acquisition speeds. However, the optical inverse problem based on such instrumental approach is still ill-conditioned. Herein, we investigate the benefit of employing compressive sensing-based preconditioning to wide-field structured illumination and detection approaches. We assess the performances of Fluorescence Molecular Tomography (FMT) when using such preconditioning methods both in silico and with experimental data. Additionally, we demonstrate that such methodology could be used to select the subset of patterns that provides optimal reconstruction performances. Lastly, we compare preconditioning data collected using a normal base that offers good experimental SNR against that directly acquired with optimal designed base. An experimental phantom study is provided to validate the proposed technique.


Introduction
Fluorescence Molecular Tomography (FMT) has known a rapid development over the last decade and half [1]. Coupled with an ever increasing availability of fluorescent molecular probes operating in the near-infrared spectra, the high sensitivity and multiplexing potential of optical systems has positioned FMT as a powerful molecular imaging technique for preclinical studies. Recently, a new compressive instrumental approach has been proposed for preclinical FMT in which wide-field structured light is employed to enable fast acquisition of tomographic data sets over large volumes. Using light modulators, illumination (and detection) bases can be created with great flexibility. Belanger et al [2] demonstrated that diffuse optical tomography based on a double light modulator architecture was achievable at very high speed. In this implementation, both illumination and detection masks are jointly employed, leading to a tomographic "single pixel" system. We successfully extended this approach to time-resolved FMT by implementing time-resolved wide-field structured light illumination and generating detection patterns via spatial integration of gated ICCD data sets. First, we demonstrated that the spatial integration did not compromise temporal data sets and that quantitative accurate optical tomography was feasible with performance similar to punctual optode strategies for absorptive inclusions [3]. Then, we demonstrated that wholebody small animal imaging FMT was feasible within a time frame of a few minutes when using such approach with a quantized low-frequency imaging base of 36 patterns [4] and cross-validation with micro-CT. Similar compressive imaging work in FMT was reported later on by Ducros et al who used a wavelet base [5]. Beside validations in well-controlled settings, the method has been applied in vivo to quantify FRET signals [6], with applications in drug delivery monitoring [7].
Recently, we have extended the methodology to hyperspectral time-resolved imaging by coupling the double light modulator design with a time-resolved spectrophotometer [8]. This implementation enables the acquisition of dense 4D data cubes for FMT (2D spatial, temporal and spectral). However, although all these compressive implementations based on wide-field structured light yield shorter acquisition time, better SNR, and less hyper-sensitivity to surface, the resulting optical inverse problem is still ill-conditioned. Sparsity constraints can be employed to improve the reconstruction distribution, especially when using early gates [9], but the inverse problem nonetheless suffers from high sensitivity to noise.
Herein, we investigate the benefits of using compressive sensing-based preconditioning to the wide-field FMT method. More precisely, we implement the preconditioning technique proposed by Jin et al [10] that was devised for point source/detection strategies to the widefield illumination/detection strategy. The merit of the technique is evaluated when applied to the whole Jacobian or independently to the illumination/detection bases. Beyond applying this preconditioning method, we propose a novel approach to identify the subset of preconditioned patterns that leads to the optimal reconstruction results. Last, we compare results when using the preconditioning method post-hoc on experimental data from our bar pattern base with that directly collected from the optimal base.

Method
Optical tomography is an ill-conditioned inverse problem that is notoriously difficult to solve [11]. In the case of FMT, where fluorescence molecular probes are employed, thankfully the imaging space is, by design, sparse. Hence, there has been an increasing interest in using compressive sensing based methods to enhance FMT performance. To date, the main focus of compressive sensing based approaches has been related to the use of sparsity constraints when solving the inverse problem iteratively [9,[11][12][13][14][15]. Though, compressive sensing methods can also be devised to optimize the Jacobian features for improved sparse signal recovery [16][17][18]. Recently, Jin et al [10] proposed a preconditioning method based on compressive sensing for punctual optode data sets and recovery of sparse fluorophore distributions. The approach was validated in a phantom experiment and improvements in FMT performances were demonstrated. Herein, we follow the same strategy and apply it to the wide field structured light instrumental paradigm. We provide below the main salient points of the formulation of the inverse problem, theoretical derivation of the preconditioners, inverse solvers and objective metrics employed.

Forward model and inverse problem
To cast the optical inverse problem, we follow the formulation described in [19]. In such formulation, a Monte Carlo (MC) method [20] is employed to derive the excitation and emission forward models provided by (simplified to the continuous case): where φ is light field and subscripts x, m stand for excitation and emission wavelengths.
∈ Ω r denotes location and i s denotes the i th source. a μ is absorption coefficient, axf ημ is fluorophore yield and x g / m g are Green's function computed via the MC approach. Then, the measurement at the j th detector due to i th source can be expressed as: A linear relationship between measurements b and fluorophore yields x can be derived as: where A is the sensitivity/Jacobian matrix, x is the discretized image vector, and b is the vector of measurements. M is the total number of source-detector pairs and N is the number of elements in the imaging domain. The matrix A is always ill-conditioned and can be further manipulated (preconditioned) for improved reconstructions.

Compressive sensing-based preconditioning
In this section, we provide a brief summary of the theoretical derivation of the preconditioners which follow the formulation described in Ref [10]. The sensitivity matrix A is actually the column-wise Kronecker product of the excitation light field Φ via structured illumination and the emission light field G acquired via wide-field detection: in which ( ) . Reducing these norms can be performed via preconditioning Φ and G where s M and d M are the wide-field structured light excitation and detection preconditioners, also known as optical and measurement masks or separate masks.
Alternatively, we can also apply a preconditioning matrix A M , also called global mask, to the whole Jacobian to reduce its coherence. Following the derivation from Jin et al [10], the preconditioning matrices can be obtained via Singular Value Decomposition (i.e. Φ = UΣV ). The mathematical expressions for these preconditioners are similar and expressed as: where i stands for the respective preconditioner s, d or A, T i i i Λ = Σ Σ and ε>0 is a regularization parameter for stabilization in case of poor condition number. Lastly, since s M and d M might have negative entries and cannot be directly applied for experiments, they are decomposed into two matrices with non-negative values such as:

Reconstruction algorithms and objective evaluation metrics
The performance of FMT is highly dependent on the selection of the inverse solver and parameters. Herein, we employ two different solvers. First, for all in silico investigations, we employ a popular solver, the Matlab built-in LSQR function. As we focus solely on assessing the merits of the CS-based preconditioning method and numerous tests are performed, this solver provides a fast and reliable approach for comparing different preconditioning strategies. The convergence tolerance is set to 10 −4 and the iteration number is set equal to the noise level. For example, if we add 25dB white Gaussian noise, the iteration number will then be 25. No regularization scheme is employed jointly with the LSQR solver. However, LSQR is not expected to provide the best reconstruction results for our application. Hence, for experimental data, we also employ our previously reported L p regularization method [9], in which L p (0 < p ≤ 1) norm is added to the optimization problem as sparsity constraint. Such solver is expected to be optimal for recovery of sparse signals. An L-curve is plotted in all cases (L 1 ) to objectively select the best regularization parameter [9,15]. L-curve analysis was performed using the same approach described in Ref [21].
To objectively compare the tomographic reconstruction quality, we established three quantitative metrics: localization error (LE), volume error (VE) and root mean square error (RMSE) under isovolume set at 50%. LE is defined as the Euclidean distance between the 3D centroid of the reconstructed object (within the isovolume) and the ground truth, in units of mm. VE is defined as the difference of reconstructed object volume (within the isovolume) and the ground truth volume divided by the volume of ground truth.
RMSE is a good metric for reporting the quantitative accuracy of the reconstructions and is defined as: where truth X is the ground truth object and recon X is the reconstructed object. N is the total number of nodes in the imaging domain.

Numerical set up
For the in silico study, we built a numerical phantom of size 50mm × 40mm × 20mm with optical properties set to μ a = 0.05 cm −1 , μ s ' = 8 cm −1 and g = 0.9 (Henyey-Greenstein phase function). Two fluorescence inclusions with the same fluorescence yield, length (20mm), and radius (3mm), were positioned in the middle of the phantom (at a height of 10mm). The two inclusions were set 15mm apart (see Fig. 1(a)). The overall imaging domain was discretized with a homogenous mesh containing 9,678 nodes.
To replicate our typical experimental settings, we simulated quantized illumination and detection pattern bases comprised of 40 bar-shape patterns with overall size of 40mmx40mm ( Fig. 1(b)). This base led to the best imaging performances as reported in [2]. Additionally, such base is easy to implement for fast acquisition and does not depend on a priori knowledge of the sample parameters. If such parameters, i.e. animal size/posture/geometry and optical properties, are known prior to the imaging session, then the optimal base can be derived using the approach proposed by Dutta et al. [22].
Measurements in transmission geometry were obtained for all combinations of illumination and detection patterns (1,600) using our wide-field mesh-based MC code [23]. This code allows efficient simulation of complex extended illumination sources with accuracy and construction of the associated Jacobians using the adjoint method [20]. 10 8 photons were used per simulated pattern. Note that the simulations were performed in the time domain and then integrated to yield continuous-wave (CW) data sets. Different levels of noise were simulated by adding Gaussian noise to the synthetic tomographic data sets. The simulated Signal-to-Noise ratios (SNR) ranged from 20dB to 40dB, which was relevant to experimental noise typically encountered in FMT.

Masks and resulting patterns
Applying the preconditioning method as described in section 2.2 yields potentially three masks: illumination, detection and global masks. In all cases, these preconditioners are computational constructs, though in the case of illumination and detection masks, they can also be applied experimentally. These masks, when applied to the quantized low-frequency imaging base, yield an optimal base. Each optimal pattern is computed as the weighted superposition of 40 original patterns, with weights provided in one row of the preconditioner. Each row of the preconditioner then leads to one optimal structured light illumination/ detection pattern. As each separate mask is decomposed into two matrices with non-negative entries, the optimal illumination and detection patterns come in pairs. Figure 2 provides a visual depiction of the paired preconditioning matrices and some associated optimal patterns. It is worthy to note that masks and resulting patterns are completely independent of the object being imaged, so no priori information is required to apply the proposed scheme.

Effect of regularization parameter on the preconditioning
As described in Eq. (8), a regularization parameter ε is implemented in the SVD-driven derivation of the preconditioner. Such a regularization parameter is necessary to obtain a stable matrix due to poor conditioning of the illumination/detection structured wide-field light field matrices, as well as the full Jacobian. The regularization parameter plays a critical role in regulating matrix orthogonality, noise amplification, and high-frequency information. To provide insights on the effect of this parameter, we provide in Fig. 3 the ranked values of the normalized inner product of the Jacobian with global preconditioning ( A M ) and separate masks ( s M , d M ). As expected, increasing the regularization parameter value will decrease the matrix orthogonality, impairing the preconditioning effect. Hence, low regularization parameters are preferred to retain information content. However, noise propagation can then become an issue. In Fig. 4, we provide a pictorial description of this effect. In the case of perfect (noiseless) measurements, the noise propagation is limited even in the case of ε = 10 −3 (Fig.  4(d)) whereas it becomes overwhelming at noise of 25dB (Fig. 4(h)). The same trend was observed in the reconstruction results (not shown) in which use of a larger regularization parameter reduced the localization error and root mean square error whereas the volume error stayed large or increased since high regularization is associated with loss of resolution. Based on these results, ε = 10 −2 was selected for subsequent studies.

Selection of pattern/measurement subset
We have seen from the last section that noise is over-amplified for some source-detector pairs and could impair reconstruction results. We also notice that, after preconditioning the information content is mainly associated with a few patterns (top left corner as shown in Fig.  4). Hence, selecting this subset of measurements after preconditioning should provide the most useful and accurate sub data set. Herein, we propose to select this sub data set in two steps. From Eq. (8) we can see that every preconditioning pattern (k th ) is associated to one singular value (k th ) in the ranked SVD. As it is well-established that most relevant measurements will be requested by the first m highest singular values, we can select the associated patterns to form a sub-base that contains the most important and relevant information. This process typically decreases the number of patterns employed from couple of dozens to around ten. In step two, measurements with the largest absolute values are picked out for reconstruction. An example of the process of pattern and measurement selection is provided in Fig. 5. The size of the sub data set depends on multiple parameters such as the experimental base used, fluorophore distribution and SNR. A few sub data sets should be tested and the one leading to the best reconstruction is selected.

Comparison of separate masks and global mask preconditioning
To establish the merits of the preconditioning technique, reconstructions of the numerical phantom with separate masks, global mask, and no preconditioning were performed at levels of noise ranging from 20dB to 40dB with 5dB steps. At each level of noise, 100 independent synthetic tomographic data sets were created and reconstructed (a total of 1,500 reconstructions). For both preconditioning strategies, 7 patterns/49 first source-detector pairs and the 13 largest measurement values are selected for reconstruction. Each individual reconstruction was post-processed to yield LE, VE and RMSE. We report in Table 1 the mean values computed from these 100 trials for each condition investigated (noise level, preconditioning used). As discussed in section 3.4, the optimal pattern/measurement number varies from case to case, so not all reconstructions for global mask/separate masks preconditioning in the 100 reconstructions provide the best results. Even so, preconditioning leads to improved reconstructions versus using no preconditioning when considering RMSE and LE in all cases. In terms of VE, the performance is a little worse due to loss of high-frequency information when selecting the pattern/measurement subset. Furthermore, the reconstructions are performed using a simple LSQR algorithm without use of regularization. We propose in Fig.  6 the visualization of a reconstruction example for 25dB noise. While both types of preconditioning improve reconstructions by reducing artifacts and providing a better rendering of the structural details, only separate mask preconditioning is able to fully resolve two objects, which explains the smaller VE for separate masks than global mask preconditioning. The difference could be due to larger condition number and thus more significant noise amplification (the condition number of the sensitivity matrix after applying separate masks is 5.4 × 10 6 while it's 3.8 × 10 7 for global mask). This is in agreement with observations reported in ref 10. Hence, we conclude that separate masks preconditioning provides better performance in terms of reconstructions than global mask preconditioning.

Comparison of preconditioning with separate masks and imaging with optimal base
If preconditioning yields better results, we can hypothesize that in conditions where phantom properties (geometrical and optical) are known, it is better-suited to compute ad hoc the optimal structured illumination and detection bases and directly use them experimentally. To do so, a positive and negative base, as shown in Fig. 2, can be constructed to sequentially acquire measurements. This approach is classical in CS-based imaging and has been implemented for wavelet bases in FMT by Ducros et al. [24] (termed virtual sources). The caveat of this approach is that it requires at least twice the acquisition time since the positive and negative components need to be acquired sequentially to form the mathematical base. When directly applying the optimal pattern base, the measurements used for reconstruction are computed from four experimental measurements, such as where m 1 is collected from i P + and d P + , m 2 from i P + and d P − , m 3 from i P − and d P + , m 4 from i P − and d P − . i P + , i P − stand for positive/negative optimal illumination patterns and d P + , d P − stand for positive/negative optimal detection patterns.
To assess the merits of directly imaging using the optimal illumination/detection base as obtained from SVD or performing preconditioning on experimental data from the quantized low-frequency base, we propose in Fig. 7 the measurements associated with each case at different noise levels. SNR reported are associated with each measurement (m i ) in the case of optimal base imaging, and for each measurement as in Fig. 4(e) in the case of quantized low frequency base. Fig. 7. Single-pixel measurements for first 81 pairs of illumination and detection patterns in the case of preconditioning after imaging with quantized low-frequency base (first row) and directly imaging with optimal base (second row). From left to right, 50dB, 40dB, 30dB and 20dB Gaussian noise is added.
As seen, directly employing the optimal base leads to very poor SNR in the measurements whereas the preconditioned measurements using quantized low frequency base experimental data are robust even at 20dB. This is expected as optimal base imaging data are computed out of 4 experimental measurements whereas preconditioned measurements using quantized low frequency base experimental data are obtained by the weighed sum of 1,600 measurements. Hence, these results indicate that it is preferable to acquire data sets with high SNR values using a quantized base and then perform preconditioning instead of using a mathematically optimal base with poor SNR. Note also that the optimal base as described in Visualization 1 are challenging to implement experimentally. Using the optimal positive illumination pattern set as an example, the range (difference between maximum and minimum pixel value) is between 100 and 150 for pattern 1-9 (out of the 0-255 potential values using the 8-bit DMDs typically employed), 60-80 for pattern 10-21 and almost 0 for the rest. Hence, it is impossible to experimentally use generated patterns after #21. In addition, for the first 9 patterns, the minimum difference between pixel values could be as low as 0.003 (3rd pattern), well below the minimum value that can be encoded using our DMDs.

Experimental validation
A 20mm thick phantom with optical properties of μ a = 0.1cm −1 and μ s ' = 9.6cm −1 was created using the recipe described in ref 6 and validated via time-resolved spectroscopy. Two tubes (3.0mm diameter and 15.5mm center-center distance) with different fluorescence dyes, Alexa Fluor 750 and 3.3′diethylthiatricarbocyanine iodide, were placed at 10mm depth in the phantom. At both illumination and detection side, light is encoded by a digital mirror device (DMD) over 40mm by 30mm area on the sample with 36 quantized low frequency patterns as used in [8]. Excitation is performed at 750nm and emission light is collected through a TCSPC module over 10 wavelength channels (centered at 770nm-820nm) and 256 time gates (48.9ps width). We integrated the data temporally and over 3 spectral channels (channel 775-785) to obtain a continuous-wave data set. We applied the preconditioning approach to this data set. Figure 8 provides the experimental data set as well as the measurements after preconditioning. The relative regularization parameter was set to ε = 10 −3 due to the high SNR level of data acquired; 7 patterns and 25 measurements after preconditioning were selected. Reconstructions were performed using both LSQR and L 1 -regularization methods, with data set before and after preconditioning. The number of iterations is set to 30 for LSQR in cases and 3 × 10 4 , 3 × 10 5 for L 1 -regularization before/after preconditioning. The regularization parameters were chosen from the L-curves. The reconstructions results are provided in Fig. 9. In all cases, separate mask preconditioning led to better reconstructions. Especially, when associated with sparsity constraints, the two inclusions were separated efficiently and located with accuracy. Note that these reconstructions were performed using a CW data type and better performances are expected when using early gates [9].

Conclusion
We have implemented and investigated a compressive sensing based preconditioning method for wide-field structured light fluorescence molecular tomography. The method was tested in silico and experimentally. In all cases investigated, separate mask preconditioning, in which structured illumination and detection light fields are preconditioned independently, led to the best overall performances in 3D reconstructions. Additionally, we demonstrated that imaging with a quantized low-frequency base, in conjunction with preconditioning, led to better results than directly imaging with optimal base as derived from the SVD.
However, our investigation was limited to low-frequency quantized basis and the CW data type. For future studies, we plan to apply this preconditioning scheme to 1) various kinds of commonly used basis in compressive imaging, such as wavelet based or DCT based patterns; 2) time-domain imaging, especially for early gates, which offer improved resolution. Also, we plan to extend the method to hyperspectral data sets. 4D data cubes lead to large data sets that are difficult to employ for FMT due to memory constraints. The selection of pattern/measurement subset approach described herein, leads to significantly reduced data sets that would allow us to leverage the rich information content captured by compressive hyperspectral time-resolved imaging for fast whole-body multiplexed FMT.