Improving mesoscopic fluorescence molecular tomography via preconditioning and regularization

Fugang Yang; Ruoyang Yao; Mehmet Saadeddin Ozturk; Denzel Faulkner; Qinglan Qu; Xavier Intes

doi:10.1364/BOE.9.002765

1. Introduction

As a novel imaging modality combining the benefits of millimeter deep molecular imaging with high sensitivity, Mesoscopic Fluorescence Molecular Tomography (MFMT) can resolve the distribution and concentration of fluorophores inside biological tissues based on boundary measurements [1,2]. MFMT has the unique potential to cover the gap between microscopic and macroscopic examination of biotissues by enabling resolution up to hundred micrometers level at a depth of a few millimeters in intact tissues (without optical windows or clearing agents). MFMT has already found applications in imaging bio-printed thick constructs [3,4], brain molecular [5,6], dental imaging [7] and tumor xenografts studies [8,9]. However, performances of MFMT are intrinsically associated with the inverse optical problem that need to be solved.

Typically, the MFMT inverse problem is ill-posed and ill-conditioned. Hence it is very challenging to solve it in a robust fashion. Among the various techniques developed to solve ill-conditioned problems, iterative methods, embedded with various optimization algorithms, have become an effective approach that use successive approximations to obtain solutions with the help of additional constraints introduced by appropriate regularization techniques [10,11]. However, the performance of the iterative method adopted depends greatly on the incoherence and the spectrum of the forward sensitivity matrix. Additionally, many iterative methodologies rely on preconditioners to improve performance and ensure fast convergence [12]. In the case of the MFMT inverse problem, the measurements acquired are characterized with high redundancy due to the diffuse nature of the light collected, reflection geometry employed and dense spatial sampling. Hence, the associated Jacobian demonstrates high coherence between its columns and high condition number due to fine discretization of the volume to be imaged to attain hundred micrometers resolution. In turn, the designed iterative optimization algorithm will slowly converge and may be susceptible to the propagation and amplification of errors after successive iterations [13].

On the other hand, Compressed Sensing (CS) methodologies have been proved to be successful methods to recover sparse signals from far fewer samples than required by the traditional Shannon-Nyquist sampling theorem [14,15]. The CS framework guarantees accurate recovery of sparse signals under certain conditions, one of which is incoherence, i.e. in our case, the degree of the orthogonality of forward sensitivity matrix. Recent works in diffuse optical tomography (DOT) have demonstrated that using a properly preconditioning matrix can reduce the coherence between columns of the forward sensitivity matrix, and as a consequence alleviate the ill-posedness of the underdetermined linear system [16,17]. Similarly, our group has proposed a preconditioning strategy to reduce the coherence of sensitivity matrix for wide-field fluorescence molecular tomography [18]. As for the MFMT inverse problem itself, Yang et al. [19] investigated the influence of noise on reconstruction performance and addressed a two-step data reduction approach to achieve high fidelity results. Nevertheless, the implementation of effective preconditioning methodologies in MFMT has not been reported to date but is expected to greatly improve the performances of the reconstruction process.

In this paper, we propose two computational methods to improve both incoherence and condition number of the sensitivity matrix of MFMT. First, we propose a preconditioner to reduce the correlation between columns of the Jacobian, which contribute to mitigating the ill-posedness of the inverse problem. Then another well-designed regularization parameter is added to improve the convergence rate of the inverse solver to the targeted tolerance. The proposed strategies are described in Section 2 in details. Section 3 presents the objective metrics to evaluate spectrum property, coherence of Jacobian, noise level of measurements, and the MFMT reconstruction performances. Section 4 summarizes the in silico and in vitro experiment setups with reconstruction results. The discussion and conclusion are provided in Section 5.

2. Method

2.1 Configuration of the MFMT optical system

The optical configuration of our second generation MFMT system has been described in details in previous reports [4,19]. Here, we provide the main salient features that are relevant to the optical inverse problem. Briefly, the system is built around the raster scanning of an illumination spot over the specimen with detection performed by an electron-multiplying CCD camera (emCCD) (iXon^EM+ DU-897 back illuminated, Andor Technology) acquiring data in a descanned configuration. This set up enables the acquisition of up to 512 × 512 measurements in parallel per illumination spot, leading to very dense spatial data sets. To obtain a data set that can be inverted efficiently, the emCCD data can be binned into super-pixels that still provide high accuracy and resolution [19]. The discretization of the image space and positioning of the optodes configuration are illustrated in Fig. 1, where the green squares, red squares and blue squares respectively represent the discretized voxels of the imaging volume, pixels of emCCD camera, and super-pixels binned at 2 × 2 to provide enhanced SNR.

Fig. 1 Schematic diagram of scanning mode of the 2^nd generation MFMT system. (a) displays the discretization of imaging volume. (b) shows the dimension of the emCCD with super pixels (blue squares) binned by 2 × 2 pixels(red squares) used as detectors. (c) and (d) illustrate the scanning trajectory (blue dash) from start scanning spot (c) to the last one (d).

Download Full Size | PDF

This configuration and on the chip binning strategy leads to a total of 256 × 256 super pixels as detectors that cover a detection FOV of 6.4 × 6.4 mm². To cast the inverse problem, we further down sample the measurements space by selecting a 7 × 7 detector matrix (D₁, D₂, ……, D₄₉) with a separation of ~0.6mm. To improve the dynamical range, the central detection which coincides with the position of the illumination spot is occluded during acquisition so detector 25 (D₂₅) is not employed. On the illumination side, the raster scanning step size and dwelling time are typically set to 100 µm, 20 ms respectively, leading to a total of 961 scanning locations in the illumination FOV of 3.1 × 3.1 mm². Last, the volume to be imaged is discretized uniformly in voxels of 100 μm length, leading to R_i (i = 1, ……, 31), C_j (j = 1, ……,31), and Z_k (k = 1, 2,……, 30) number of voxels along the x, y and z axes respectively. Overall, this configuration leads to 46,128 measurements that are employed in the inverse optical problem (out of 251,920,384 possible without binning).

2.2 Forward model and inverse problem

The source-detector separation in our MFMT implementation can be as small as 600 μm. Hence, in view of typical optical properties of tissues, the collected photons cannot be accurately modeled by the Diffusion Equation [20]. Therefore, we employ Monte Carlo (MC) method, commonly known as the gold standard [21]. The MC method is well established as being accurate in modeling photon propagation in shallow tissues and near sources [22,23]. To meet the accuracy requirements under high resolution of MFMT operating in the mesoscopic regime, up to 10⁸ photons are simulated for each source/detector with the GPU-based software MCX to simulate continuous-wave (CW) Green’s functions G^x and G^m at excitation and emission wavelengths [24]. The weight function W for one source-detector pair is then computed with an adjoint formulation for efficiency [25]:

W(r_{s}, r_{d},r)= G^{x} (r_{s},r) \times G^{m} (r, r_{d})

where r_s, r_d, and r are locations for source, detector and any location inside the specimen being probed, respectively. The measured fluorescence intensity for the corresponding source-detector pair U(r_s, r_d) can be formulated as the integral equation below:

U(r_{s}, r_{d})= \int_{Ω} W(r_{s}, r_{d},r)η(r)d r^{3}

where η(r) is the 3-D distribution of the fluorophore’s effective quantum yield. As at this time our experimental system is limited to the continuous wave (CW) data type for MFMT applications, however it’s straightforward to expand both Eqs. (1,2) to time-dependent form given time-resolved Green’s functions and the lifetime of fluorophore τ [26]. With M source-detector pairs, and N discretized voxels of the Region of Interest (ROI), the inverse problem for MFMT in the matrix form becomes:

AX=b

where A∈R^{M × N} is the Jacobian, X∈R^{N × 1} represents the effective quantum yield, b∈R^{M × 1} is the vector of detector readings at all source positions.

The CS based framework has been demonstrated as an effective methodology to exactly retrieve sparse solutions from under-determined linear systems [10,11,14]. Fortunately, the solutions in our MFMT application are inherently sparse because fluorescent biomarkers are supposed to be clustered within small regions of tissues such as tumors. To achieve sparsity in the image space, we add the commonly used l₁-norm regularization term so that the optimization problem now becomes:

min {{‖ AX-b ‖}_{2}^{2} +λ {‖ X ‖}_{1}

where

{‖ X ‖}_{1} = \sum_{1}^{n} | x_{i} |

denotes the l₁-norm of fluorophore distribution vector and λ is the regularization parameter. The minimization problem in Eq. (4) can then be transformed to a convex quadratic problem with linear inequality constrains, and solved by a variety of inverse solvers such as the interior-point method that we applied in this work [27].

2.3 Preconditioning to reduce the coherence of the sensitivity matrix

The first preconditioning strategy applied aims at reducing the coherence of the sensitivity matrix, so that the performance of CS-based sparse signal recovery algorithms is improved. The design of preconditioner follows a similar approach as described in [16–18] so we only briefly introduce the derivation. When a preconditioner M_A is applied, the original Jacobian matrix becomes:

A_{pre} = M_{A} A

For simplicity, we assume that the preconditioned sensitivity matrix A_pre is also column-normalized, and to minimize the coherence of A_pre, we seek to determine M_A such that the Gramian matrix $A_{pre}^{T} A_{pre}$ approximates the identity matrix [10,11,28],

A_{pre}^{T} A_{pre} = A^{T} M_{A}^{T} M_{A} A \approx I

Equation (6) still holds after multiplying A on the left and A^T on the right,

A A_{p r e}^{T} A_{p r e} A^{T} = A A^{T} M_{A}^{T} M_{A} A A^{T} \approx A A^{T}

Taking the Singular Value Decomposition of A = U_A∑_{A $V_{A}^{T}$}, and substitute it into Eq. (7), we get:

(Σ_{A} Σ_{A}^{T}) U_{A}^{T} M_{A}^{T} M_{A} U_{A} (Σ_{A} Σ_{A}^{T}) \approx Σ_{A} Σ_{A}^{T}

Thus, we choose M_A as the preconditioner,

M_{A} = {(Σ_{A} Σ_{A}^{T})}^{- \frac{1}{2}} U_{A}^{T}

During the actual implementation, a regularization parameter is usually added for stabilization when having a large condition number, so the final form of our first preconditioner becomes:

M_{A} = {(Σ_{A} Σ_{A}^{T} + ϵ I)}^{- \frac{1}{2}} U_{A}^{T}

where ϵ = 0.1 is empirically selected to provide the best trade-off between reducing matrix coherence and suppressing noise propagation in the measurements under the experimental noise level (SNR = 4) [18].

Note that although the construction of MFMT Jacobian also follows the adjoint form, we cannot separately reduce the coherence of excitation and emission matrices and compute the preconditioned sensitivity matrix through the Kronecker product, referred as “separate masks” in [18]. This is because the excitation and emission light fields of MFMT are coupled with each other, i.e., the source and CCD camera always move together and remain the same relative position as the scanning proceeds. Thus, only “global mask” strategy can be applied to MFMT, although it’s shown that “separate masks” are slightly superior in suppressing noise and improving resolution in the reconstruction results [16,18]. Though, we can harness the symmetry of the descanned acquisition configuration that leads to a block configuration in the forward model (one source-48 detectors).

Herein, we propose a “sub-preconditioning” approach in which sub-matrices are preconditioned to overall lead to a reduction of the coherence of the whole matrix. In our case, 48 measurements are collected from each of the 961 scanning positions, so we can either perform preconditioning on 961 sub-matrices with 48 rows, or on 48 sub-matrices with 961 rows. In practice, we chose the latter one because it’s more effective at reducing the overall coherence of the Jacobian. Thus, each sub-sensitivity matrix A_p (p = 1, ……, 48) corresponds to measurements from the same detector at all scanning positions. In the following studies, preconditioning strategies will be performed both globally and block by block, referred as “whole-preconditioning” and “sub-preconditioning”, and the performance of the two approaches are compared and analyzed.

2.4 WIIM to speed-up the convergence rate of the inverse solver

The convergence rate of an iterative optimal algorithm to solve a large linear system depends mainly on the condition number of its coefficient matrix [13,29]. Unfortunately the condition number in MFMT application is usually very large due to the high correlation between different source-detector pairs resulting from dense sampling and fine discretization. Hence an effective method to reduce the condition number of the sensitivity matrix is highly desired.

To this end, here we employ a regularization approach, Weighted Iterative Improvement Method (WIIM) [30], to effectively improve the ill-conditioned situation and thus contribute to a faster convergence rate of the iterative algorithm. We first construct the iterative process as following,

(Λ + γ I) x_{k + 1} = y + γ x_{k}

where Λ = A^TA is a square matrix, y = A^Tb, I is the identity matrix, and k = 1, 2, 3, … is the iteration number. The nonzero constant γ determines the degree of improvement for the condition number of Λ. A small γ leads to minor improvement of convergence rate while a very large γ could result in stagnation of iteration or even divergence. It’s usually determined from the following empirical formula [30]:

\frac{γ}{γ + δ^{2}} \leq 0.1, (γ > 0)

where

δ

is the lower bound of all nonzero elements in the sensitivity matrix. The initial value of γ can be given according to Eq. (13)

{\begin{matrix} γ = β 10^{0.5 | \log (β) | + 1} \\ β = \min (| e i g (Λ) | \end{matrix}

where β is the absolute value of the smallest eigenvalue of matrix

Λ

and log is logarithm at the base 10. Let x_k₊₁ = x_k + e_k and r_k = y - Λx_k, we can substitute them into Eq. (11) and get:

(Λ + γ I) e_{k} = r_{k}

We can apply the same CS-based solver as we solve Eq. (3) because e_k is the difference of fluorophore concentration between two iterations and thus is also sparse. With an initial guess

x_{0}

, we can then repeatedly solve Eq. (14) and update the values of e_k, r_k, and x_k. If the last residual norm ||r^k|| reaches the targeted tolerance, or the iteration number reaches the maximum, the approximated solution is reached. Note that the iteration steps described above can be easily implemented to be embedded in any selected inverse solver. Last, a positive constraints was applied during the iterative process.

3. Evaluation metrics for an objective assessment of the proposed methods

3.1 The effect of the condition number on the precision of the solution

The condition number of the sensitivity matrix determines how much a small disturbance in the measurement b could affect the solution of the inverse problem x in Eq. (3):

\frac{‖ δ x ‖}{‖ x ‖} \leq c o n d (A) \frac{‖ δ b ‖}{‖ b ‖}

where cond(A) = ||A⁻¹|| ||A|| is the condition number. A smaller condition number could lead to results that are accurate and stable against noise, because the error magnification during iterations is minimized. Therefore, it’s beneficial to apply preconditioning methods to reduce the large condition number of the MFMT sensitivity matrix.

3.2 Coherence of the sensitivity matrix

One of the most commonly used metric in the compressed sensing literature to assess the incoherence of a matrix $A$ is the mutual coherence [16], defined as the largest normalized inner product of two different columns a_p and a_q:

M (A) = \overset{m a x}{p, q, p \neq q} \frac{| 〈 a_{p}, a_{q} 〉 |}{{‖ a_{p} ‖}_{2} {‖ a_{q} ‖}_{2}}

M(A) will be small only if all the columns of A are almost orthogonal to each other, so it reports on the worst orthogonality of the whole matrix. Thus, to provide a more robust assessment of the incoherence of the sensitivity matrix, the cumulative coherence or the Babel function can be applied to report on the “average” coherence instead [31].

M_{1} (k, A) = \overset{m a x}{p} \overset{m a x}{| Q | = k, p \notin Q} \sum_{p \in Q} \frac{| 〈 a_{p}, a_{q} 〉 |}{{‖ a_{p} ‖}_{2} {‖ a_{q} ‖}_{2}}

Where |Q| is the cardinality, i.e. the number of columns, of the selected subset of the sensitivity matrix Q, and M₁(k, A) is monotonically non-decreasing [31].

3.3 Signal to Noise ratio

Noise is an important factor impacting the accuracy of the MFMT inverse problem. To describe the noise level of the MFMT system more faithfully, we define the Signal-to-Noise Ratio (SNR) in our experiment as follows:

SNR = \frac{μ (| S_{f} - S_{b} |)}{σ (| S_{b} - S_{r} |)}

where S_f is the fluorescence signal from all source positions, and S_b is the average of S_f with the ROI. S_r is also the fluorescence signal but with the sample replaced by a beam dump. So

μ (| S_{f} - S_{b} |)

indicate the mean of the true signal intensity, and

σ (| S_{b} - S_{r} |)

indicate the standard deviation of the noise due to sample readout.

For in silico experiments, a Gaussian noise with mean of zero and standard deviation of SNR is added to the simulated data set to mimic experimental measurements. We use SNR values ranging from 2 to 6 to explore the robustness of the proposed computational approaches on different noise levels.

3.4 Quantitative metrics to assess 3D reconstructions

For consistency, the same normalized figures of merit as [19] are used to quantify the difference between reconstruction results and the ground truth. Four normalized similarity metrics, normalized sum squared difference (nSSD), normalized sum absolute difference (nSAD), normalized disparity (nD), and normalized correlation (nR), are defined as follows.

nSSD (A, B) = 1 - \frac{1}{N} \sum_{i = 1}^{R} \sum_{j = 1}^{C} \sum_{k = 1}^{Z} {[A (i, j, k) - B (i, j, k)]}^{2}

nSAD (A, B) = 1 - \frac{1}{N} \sum_{i = 1}^{R} \sum_{j = 1}^{C} \sum_{k = 1}^{Z} | A (i, j, k) - B (i, j, k) |

nR= \frac{\sum_{i = 1}^{R} \sum_{j = 1}^{C} \sum_{k = 1}^{Z} A (i, j, k) \times B (i, j, k)}{\sqrt{\sum_{i = 1}^{R} \sum_{j = 1}^{C} \sum_{k = 1}^{Z} A {(i, j, k)}^{2} \times \sum_{i = 1}^{R} \sum_{j = 1}^{C} \sum_{k = 1}^{Z} B {(i, j, k)}^{2}}}

nD = 1 - \frac{1}{N} \sum_{i = 1}^{R} \sum_{j = 1}^{C} \sum_{k = 1}^{Z} A^{'} (i, j, k) \oplus B^{'} (i, j, k)

where A(i, j, k), B(i, j, k) are the value of the normalized numerical model and reconstructed result at coordinates (i, j, k), N = R × C × Z is the total voxels in a cubic phantom and R, C, Z are the number of voxels three axes, respectively.

A^{'}, B^{'}

are binaries (0 or 1) of A, B through thresholding, and

\oplus

represents the exclusive or (XOR) operation. The values of these four metrics are between 0 (for two absolute different models) and 1(for two coincident models). In other words, a larger value of the metrics indicates higher similarity between the reconstruction result and the ground truth.

4. In silico and experimental validation

4.1 Simulation settings

A numerical phantom was designed to evaluate the performance of the proposed postprocessing methods, as shown in Fig. 2. The imaging domain had a surface area of 3.1 × 3.1 mm² with a depth of 3 mm and was uniformly discretized into 31 × 31 × 30 voxels with 100 μm resolution. A vascular tree with a main trunk and three groups of offshoots was placed within the phantom at z = 1 mm. The diameter of the trunk and offshoots were 400 μm and 200 μm separately, and the separation distance between two adjacent off shoots was 100 μm. The fluorophore concentration was assumed to be uniform in the vessel with effective quantum yield equal to 1. The optical properties of domain were assumed to be homogeneous at the excitation wavelength, with absorption coefficient μ_a = 0.02 mm⁻¹, scattering coefficient μ_s' = 1 mm⁻¹, index of refraction n = 1.34, and anisotropy factor g = 0.81. These values are derived from the collagen scaffold typically employed in our bio-printing application at 6-9 mg/ml density [32] and μ_s' is also on the same level with many biological tissues [33].

Fig. 2 The numerical phantom designed to mimic a bio-printed vascular channel and sprouting capillaries. (a), (b) and (c) are the full view, xy view, and xz view of the phantom, respectively.

Download Full Size | PDF

We replicated the imaging system configuration as used in a real experiment: 961 scanning positions and 48 detectors at each source location. The in silico measurements are generated by multiplying the sensitivity matrix with the bio-printed vasculature model. We added Gaussian noise with SNR as low as 2. Then reconstructions are performed as described in Section 2. For each of the following simulations and experiment cases, the optimal regularization parameter is chosen through L-curve analysis [34].

4.2 Coherence reduction via preconditioning

We first evaluated the performance of the preconditioning method on improving the incoherence of the sensitivity matrix and reconstruction results with the numerical phantom. As we mentioned in Section 2.3, it could be directly applied on the whole sensitivity matrix (whole-pre) or block by block (sub-pre). We compared the distribution of largest 30% of normalized inner products between two columns of the sensitivity matrix in Fig. 3(a), and the cumulative coherence of the sensitivity matrix as a function of column number k in Fig. 3(b). As shown in the insert in Fig. 3(a), the sub-pre reduces the relative Area Under the Curve (AUC) to 45.27% while whole-pre reduces the AUC to 34.79%. The average slope listed in the insert in Fig. 3(b), a quantitative index of cumulative coherence, also shows a slower increment after applying the preconditioner on whole sensitivity matrix and sub sensitivity matrix, which demonstrates both preconditioning could effectively improve the orthogonality of the matrix.

Fig. 3 Curves of the normalized products (a) and cumulative coherence (b) before and after applying the preconditioner on whole sensitivity matrix and block by block sensitivity matrix.

Download Full Size | PDF

We then compared the reconstruction results between non-pre, whole-pre and sub-pre, at different measurement noise levels (from 2 to 6). The comparison of four reconstruction metrics between non-pre and whole-pre is shown in Fig. 4(a), and those between whole-pre and sub-pre is shown in Fig. 4(b). At all noise levels, the reconstruction quality with the proposed preconditioning methods outperforms that without preconditioning. Moreover, though less effective in reducing the coherence of sensitivity matrix, sub-pre can retrieve the in silico model with greater similarity and smaller disparity compared to whole-pre. The 3D reconstruction images of three cases are shown in Fig. 4(c)–4(e). In accordance with the evaluation metrics, sub-preconditioning provides visually the best result. Although whole preconditioning is better at reducing matrix coherence, the reconstruction results are less ideal than sub preconditioning. This observation is in agreement with the results in [18] and [16], which can be explained by the higher condition number and thus more severe noise amplification when applying whole-preconditioning.

Fig. 4 Reconstruction results and evaluation metrics before and after applying preconditioning. (a) shows the metrics of reconstructions of non-pre and whole-pre at different SNR levels. (b) shows the metrics of reconstructions of whole-pre and sub-pre at different SNR levels. (c-e) display the 3-D reconstruction of non-pre, whole-pre and sub-pre, respectively.

Download Full Size | PDF

4.3 Faster convergence via the weighted iteration improvement method

We further evaluated the performance of the proposed Weighted Iteration Improvement Method (WIIM) on accelerating the convergence rate of the inverse problem in four scenarios, named whole-pre without WIIM (M1), sub-pre without WIIM (M2), whole-pre with WIIM (M3) and sub-pre with WIIM (M4). It can be seen from Figs. 5(a)–5(c) that the convergence rates without WIIM always stagnate before the desired tolerance, especially for the cases with low SNR level. However, after our l₁-norm regularization based iterative algorithm is embedded with WIIM, an effective improvement on the convergence rate can be seen. The last residual norm ||r^k|| always becomes much smaller, as reported in Table 1. By comparing M1 with M3 (both whole-preconditioning) and M2 with M4 (both sub-preconditioning), we can see from Table 1 that WIIM could help improve reconstruction results for both Jacobians. On the other hand, the evaluation metrics of M3 worse than M2 shows that preconditioning has greater impact on the results than WIIM because the latter one is just a minor improvement for the iterative optimization algorithm. In addition, sub-pre converges faster than whole-pre when WIIM is applied, as illustrated by the green and blue curves in Fig. 5(a)–(c). Overall, the proposed WIIM works well for both whole sensitivity matrix and sub sensitivity matrix to fasten convergence rate.

Fig. 5 Comparison of convergence rate and visual reconstruction results among four scenarios without and with WIIM embedded in them. (a-c) plots the convergence rate and iteration number of the scenarios at three noise levels, respectively. (d-f) visualize the reconstruction of three scenarios at SNR of 4, respectively.

Download Full Size | PDF

Table 1. Quantification results of the preconditioning for numerical data (Targeted tolerance = 10⁻³, Preset max Iteration = 9,000)

View Table | View all tables in this article

The visual reconstruction results in Fig. 5(d)–(f) further demonstrate the improvement of reconstruction quality with WIIM. The obvious improvement from Fig. 5(d)–5(f) indicates that the deterioration of reconstruction could likely results from the stagnation before the iterative algorithm converges to the desired tolerance.

4.4 Performance of preconditioning methods on experimental data

Next we evaluated the performance of the proposed preconditioning and regularization methods on an experimental data set. The collagen phantom to be explored is homogeneous with optical properties μ_a = 0.02 mm⁻¹, μ_s' = 1 mm⁻¹, n = 1.34, and g = 0.81. It has a size of 3.1 × 3.1 × 3.0 mm³, and voxels of 100 × 100 × 100 μm³, same as the simulated case. Four polystyrene fluorophore beads (GFP 488/509, Cospheric) were placed 1.7~1.8 mm beneath the sample surface under standard protocols. The transversal slices across the x-y and y-z planes taken from a micro-MRI are shown in Figs. 6(a)–(b) and the whole 3-D view of the phantom, overlaid with the best reconstruction result, is shown in Fig. 6(c).

Fig. 6 Experimental phantom (a-c) are the x-y, x-z cross-sectional, and 3D view of phantom obtained from micro-MRI, overlaid with the best reconstruction result, respectively.

Download Full Size | PDF

The measurements were recorded as described in Section 2.1 and the data acquisition time at each scanning position was less than 20 ms. Similar to the in silico cases, we compared the condition number, convergence rate, and reconstruction metrics, when no / whole / sub-preconditioning is preformed, with and without WIIM. All the inverse problems were solved with our l₁-norm reconstruction algorithm as described in Section 2.2 and the optimal regularization parameters were determined through L-curves. As an example, Fig. 7(a) shows the L-curve for experimental data to determine the optimal regularization parameter when applying sub-preconditioning.

Fig. 7 Reconstruction results of experimental phantom. (a) L-curve for experimental data to determine the optimal regularization parameter when applying sub-preconditioning. (b-d) are results without preconditioning, with whole-pre and WIIM, and with sub-pre and WIIM, respectively. (e) is the result when applying sub-pre, WIIM and noise suppression on the experimental phantom .

Download Full Size | PDF

As shown in Table 2, the proposed preconditioner can effectively reduce the condition number as well as the coherence between columns of the sensitivity matrix in both formations of. It leads to a reduction of 2,524 and 3,628 times of condition number in the two sensitivity matrix formations, respectively. Meanwhile, the last residual norm also drops from 1.1479 to 0.0081with the help of WIIM, which is very close to the targeted tolerance. The four evaluation metrics of reconstruction results in Table 2 further indicate the proposed computational methods greatly improve the reconstruction fidelity. Although whole-pre and sub-pre strategies have the same last residual norm, sub-pre achieves lower condition number as well as faster convergence rate. Visual results in Fig. 7(b)–7(e) and evaluation metrics in Table 2 also validate the superiority of reconstruction quality from sub-pre. However, as shown in Figs. 7(d), even the best results still suffer from some artifact, which illustrates that the proposed methods may still be sensitive to noise. Figure 7(e) shows the reconstruction result when applying both methods on the sub sensitivity matrices and noise suppression strategy as described in [19], which is very close to the ground truth obtained from micro-MRI.

Table 2. Quantification results of the preconditioning for experimental data (Targeted tolerance = 10⁻³, Preset max Iteration = 9,000, SNR = 4)

View Table | View all tables in this article

5. Conclusion

We have proposed two computational methods to reduce the coherence and condition number of the sensitivity matrix of MFMT, and to accelerate convergence rate of the adopted iterative algorithm when solving the inverse problem. The preconditioning method can be directly applied on the whole sensitivity matrix or on sub-matrices of each detector, while the sub- preconditioning strategy provides fast convergence rate and better reconstruction quality than whole-preconditioning despite larger coherence between columns of the sensitivity matrix. This is in agreement with the observations in [16–18], because direct preconditioning usually lead to a large condition number and amplification of noise in the measurement vector. The regularization method based on WIIM proves to be effective to avoid stagnation before the iterative algorithm reaches the target tolerance. We have tested the performance of these two methods on in silico as well as in vitro data sets. In both cases, the reconstruction fidelity is significantly improved compared to previous results. We plan to further investigate the potential of the proposed methods for augmented data sets, such as in time-resolved cases [26] when applied to FRET tomography [35] to enable to monitor cellular processes [36,37] and in the case of phased array systems [38–41].

Funding

National Institutes of Health (NIH) (R01-EB019443 and R01 BRG-CA20772); Natural Science Foundation of Shandong Province (ZR2018MF034); National Natural Science Foundation of China (61472227 and 81771532).

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References and links

1. E. M. C. Hillman, D. A. Boas, A. M. Dale, and A. K. Dunn, “Laminar optical tomography: demonstration of millimeter-scale depth-resolved imaging in turbid media,” Opt. Lett. 29(14), 1650–1652 (2004). [CrossRef] [PubMed]

2. M. S. Ozturk, V. K. Lee, L. Zhao, G. Dai, and X. Intes, “Mesoscopic fluorescence molecular tomography of reporter genes in bioprinted thick tissue,” J. Biomed. Opt. 18(10), 100501 (2013). [CrossRef] [PubMed]

3. M. S. Ozturk, C. W. Chen, R. Ji, L. Zhao, B. N. Nguyen, J. P. Fisher, Y. Chen, and X. Intes, “Mesoscopic Fluorescence Molecular Tomography for Evaluating Engineered Tissues,” Ann. Biomed. Eng. 44(3), 667–679 (2016). [CrossRef] [PubMed]

4. M. S. Ozturk, X. Intes, G. dai, and V. K. Lee, “Longitudinal Volumetric Assessment of Glioblastoma Brain Tumor in 3D Bio-Printed Environment by Mesoscopic Fluorescence Molecular Tomography,” in Biomedical Optics 2016, OSA Technical Digest (online) (Optical Society of America, 2016), paper JM3A.46.

5. Q. Tang, J. Lin, V. Tsytsarev, R. S. Erzurumlu, Y. Liu, and Y. Chen, “Review of mesoscopic optical tomography for depth-resolved imaging of hemodynamic changes and neural activities,” Neurophotonics 4(1), 011009 (2016). [CrossRef] [PubMed]

6. M. Azimipour, M. Sheikhzadeh, R. Baumgartner, P. K. Cullen, F. J. Helmstetter, W. J. Chang, and R. Pashaie, “Fluorescence laminar optical tomography for brain imaging: system implementation and performance evaluation,” J. Biomed. Opt. 22(1), 16003 (2017). [CrossRef] [PubMed]

7. F. Long and X. Intes, “Dental optical tomography with upconversion nanoparticles—a feasibility study,” J. Biomed. Opt. 22(6), 66001 (2017). [CrossRef] [PubMed]

8. M. S. Ozturk, D. Rohrbach, U. Sunar, and X. Intes, “Mesoscopic fluorescence tomography of a photosensitizer (HPPH) 3D biodistribution in skin cancer,” Acad. Radiol. 21(2), 271–280 (2014). [CrossRef] [PubMed]

9. Y. Chen, S. Yuan, J. Wierwille, R. Naphas, Q. Li, T. R. Blackwell, T. W. Paul Jr, and G. Kristine, “Integrated optical coherence tomography (OCT) and fluorescence laminar optical tomography (FLOT),” IEEE J. Sel. Top. Quant. 16(4), 755–766 (2010). [CrossRef]

10. D. L. Donoho and M. Elad, “Optimally sparse representation in general (nonorthogonal) dictionaries via l minimization,” Proc. Natl. Acad. Sci. U.S.A. 100(5), 2197–2202 (2003). [CrossRef] [PubMed]

11. D. L. Donoho, Y. Tsaig, I. Drori, and J. L. Starck, “Sparse solution of underdetermined systems of linear equations by stagewise orthogonal matching pursuit,” IEEE Trans. Inf. Theory 58(2), 1094–1121 (2012). [CrossRef]

12. M. T. Schaub, M. Trefois, P. van Dooren, and J.-C. Delvenne, “Sparse matrix factorizations for fast linear solvers with application to Laplacian systems,” SIAM J. Matrix Anal. Appl. 38(2), 505–529 (2017). [CrossRef]

13. M. Benzi, “Preconditioning techniques for large linear systems: a survey,” J. Comput. Phys. 182(2), 418–477 (2002). [CrossRef]

14. F. Yang, M. S. Ozturk, L. Zhao, W. Cong, G. Wang, and X. Intes, “High-resolution mesoscopic fluorescence molecular tomography based on compressive sensing,” IEEE Trans. Biomed. Eng. 62(1), 248–255 (2015). [CrossRef] [PubMed]

15. Q. Pian, R. Yao, N. Sinsuebphon, and X. Intes, “Compressive hyperspectral time-resolved wide-field fluorescence lifetime imaging,” Nat. Photonics 11(7), 411–414 (2017). [CrossRef] [PubMed]

16. A. Jin, B. Yazici, and V. Ntziachristos, “Light illumination and detection patterns for fluorescence diffuse optical tomography based on compressive sensing,” IEEE Trans. Image Process. 23(6), 2609–2624 (2014). [CrossRef] [PubMed]

17. A. Jin, B. Yazici, A. Ale, and V. Ntziachristos, “Preconditioning of the fluorescence diffuse optical tomography sensing matrix based on compressive sensing,” Opt. Lett. 37(20), 4326–4328 (2012). [CrossRef] [PubMed]

18. R. Yao, Q. Pian, and X. Intes, “Wide-field fluorescence molecular tomography with compressive sensing based preconditioning,” Biomed. Opt. Express 6(12), 4887–4898 (2015). [CrossRef] [PubMed]

19. F. Yang, M. S. Ozturk, R. Yao, and X. Intes, “Improving mesoscopic fluorescence molecular tomography through data reduction,” Biomed. Opt. Express 8(8), 3868–3881 (2017). [CrossRef] [PubMed]

20. V. L. Wang, and I. W. Hsin Biomedical Optics: Principles and Imaging (John Wiley & Sons, 2012).

21. C. Zhu and Q. Liu, “Review of Monte Carlo modeling of light transport in tissues,” J. Biomed. Opt. 18(5), 50902 (2013). [CrossRef] [PubMed]

22. L. Wang and S. L. Jacques, “Hybrid model of Monte Carlo simulation and diffusion theory for light reflectance by turbid media,” J. Opt. Soc. Am. A 10(8), 1746–1752 (1993). [CrossRef] [PubMed]

23. S. H. Tseng, C. Hayakawa, J. Spanier, and A. J. Durkin, “Investigation of a probe design for facilitating the uses of the standard photon diffusion equation at short source-detector separations: Monte Carlo simulations,” J. Biomed. Opt. 14(5), 054043 (2009). [CrossRef] [PubMed]

24. Q. Fang and D. A. Boas, “Monte Carlo simulation of photon migration in 3D turbid media accelerated by graphics processing units,” Opt. Express 17(22), 20178–20190 (2009). [CrossRef] [PubMed]

25. J. Chen and X. Intes, “Comparison of Monte Carlo methods for fluorescence molecular tomography-computational efficiency,” Med. Phys. 38(10), 5788–5798 (2011). [CrossRef] [PubMed]

26. J. Chen, V. Venugopal, and X. Intes, “Monte Carlo based method for fluorescence tomographic imaging with lifetime multiplexing using time gates,” Biomed. Opt. Express 2(4), 871–886 (2011). [CrossRef] [PubMed]

27. S.-J. Kim, K. Koh, M. Lustig, S. Boyd, and D. Gorinevsky, “An interior-point method for large-scale l1-regularized least squares,” IEEE J. Sel. Top. Signal Process. 1(4), 606–617 (2007). [CrossRef]

28. J. M. Duarte-Carvajalino and G. Sapiro, “Learning to sense sparse signals: simultaneous sensing matrix and sparsifying dictionary optimization,” IEEE Trans. Image Process. 18(7), 1395–1408 (2009). [CrossRef] [PubMed]

29. Y. Saad and H. A. van der Vorst, “Iterative solution of linear systems in the 20th century,” J. Comput. Appl. Math. 123(1–2), 1–33 (2000). [CrossRef]

30. X. Deng, L. Yin, S. Peng, and M. Ding, “An iterative algorithm for solving ill-conditioned linear least squares problems,” Geod. Geodyn. 6(6), 453–459 (2015). [CrossRef]

31. J. A. Tropp, “Greed is Good: Algorithmic Results for Sparse Approximation,” IEEE Trans. Inf. Theory 50(10), 2231–2242 (2004). [CrossRef]

32. L. Zhao, V. K. Lee, S.-S. Yoo, G. Dai, and X. Intes, “The integration of 3-D cell printing and mesoscopic fluorescence molecular tomography of vascular constructs within thick hydrogel scaffolds,” Biomaterials 33(21), 5325–5332 (2012). [CrossRef] [PubMed]

33. S. L. Jacques, “Optical properties of biological tissues: a review,” Phys. Med. Biol. 58(11), R37–R61 (2013). [CrossRef] [PubMed]

34. L. Zhao, H. Yang, W. Cong, G. Wang, and X. Intes, “Lp regularization for early gate fluorescence molecular tomography,” Opt. Lett. 39(14), 4156–4159 (2014). [CrossRef] [PubMed]

35. V. Venugopal, J. Chen, M. Barroso, and X. Intes, “Quantitative tomographic imaging of intermolecular FRET in small animals,” Biomed. Opt. Express 3(12), 3161–3175 (2012). [CrossRef] [PubMed]

36. K. Abe, L. Zhao, A. Periasamy, X. Intes, and M. Barroso, “Non-invasive in vivo imaging of near infrared-labeled transferrin in breast cancer cells and tumors using fluorescence lifetime FRET,” PLoS One 8(11), e80269 (2013). [CrossRef] [PubMed]

37. S. Rajoria, L. Zhao, X. Intes, and M. Barroso, “FLIM-FRET for cancer application,” Curr. Mol. Imaging 3(2), 144–161 (2015). [CrossRef] [PubMed]

38. X. Intes, B. Chance, M. Holboke, and A. Yodh, “Interfering diffusive photon-density waves with an absorbing-fluorescent inhomogeneity,” Opt. Express 8(3), 223–231 (2001). [CrossRef] [PubMed]

39. Y. Chen, C. Mu, X. Intes, D. Blessington, and B. Chance, “Near-infrared phase cancellation instrument for fast and accurate localization of fluorescent heterogeneity,” Rev. Sci. Instrum. 74(7), 3466–3473 (2003). [CrossRef]

40. X. Intes, V. Ntziachristos, and B. Chance, “Analytical model for dual-interfering sources Diffuse Optical Tomography,” Opt. Express 10(1), 2–14 (2002). [CrossRef] [PubMed]

41. V. A. Markel, V. Mital, and J. C. Schotland, “Inverse problem in optical diffusion tomography. III. Inversion formulas and singular-value decomposition,” J. Opt. Soc. Am. A 20(5), 890–902 (2003). [CrossRef] [PubMed]

SNR	Method	Condition number	Residual norm	Number of Iteration	nR	nD	nSAD	nSSD
2	whole-pre without WIIM (M1)	256.42	1.2207	4,123	0.8654	0.6798	0.5245	0.4880
	sub-pre without WIIM (M2)	178.29	0.9656	4,672	0.9027	0.8091	0.5488	0.8223
	whole-pre with WIIM (M3)	42.76	0.1017	8,979	0.8727	0.6910	0.5476	0.4920
	sub-pre with WIIM (M4)	25.19	0.1017	7,226	0.9129	0.8221	0.5686	0.8316
4	whole-pre without WIIM (M1)	256.42	1.1479	3,963	0.9143	0.7822	0.6134	0.7908
	sub-pre without WIIM (M2)	178.29	0.8654	4,419	0.9327	0.8398	0.6577	0.8674
	whole-pre with WIIM (M3)	42.76	0.0081	8,085	0.9278	0.8091	0.6930	0.8242
	sub-pre with WIIM (M4)	25.19	0.0081	5,853	0.9429	0.8510	0.7123	0.8821
6	whole-pre without WIIM (M1)	256.42	0.0028	4,267	0.9443	0.8901	0.6576	0.8951
	sub-pre without WIIM (M2)	178.29	0.0026	4,978	0.9443	0.9218	0.7005	0.9241
	whole-pre with WIIM (M3)	42.76	0.00089	7,774	0.9528	0.9091	0.7212	0.9258
	sub-pre with WIIM (M4)	25.19	0.00089	4,923	0.9619	0.9315	0.7406	0.9367

Method	Condition number	Relative AUC	Residual norm	Number of Iteration	nR	nD	nSAD	nSSD
non-preconditioning	648,496	100%	1.1479	3,963	0.5440	0.6517	0.4562	0.4809
whole-pre without WIIM	268.91	34.79%	1.1426	7,910	0.9079	0.7762	0.6106	0.7902
sub-pre without WIIM	189.75	45.27%	1.1437	5,792	0.9270	0.8280	0.6513	0.8593
whole-pre with WIIM	49.08	34.79%	0.0081	8,085	0.9173	0.8075	0.6898	0.8217
sub-pre with WIIM	26.85	45.27%	0.0081	5,853	0.9379	0.8478	0.7065	0.8712

SNR	Method	Condition number	Residual norm	Number of Iteration	nR	nD	nSAD	nSSD
2	whole-pre without WIIM (M1)	256.42	1.2207	4,123	0.8654	0.6798	0.5245	0.4880
	sub-pre without WIIM (M2)	178.29	0.9656	4,672	0.9027	0.8091	0.5488	0.8223
	whole-pre with WIIM (M3)	42.76	0.1017	8,979	0.8727	0.6910	0.5476	0.4920
	sub-pre with WIIM (M4)	25.19	0.1017	7,226	0.9129	0.8221	0.5686	0.8316
4	whole-pre without WIIM (M1)	256.42	1.1479	3,963	0.9143	0.7822	0.6134	0.7908
	sub-pre without WIIM (M2)	178.29	0.8654	4,419	0.9327	0.8398	0.6577	0.8674
	whole-pre with WIIM (M3)	42.76	0.0081	8,085	0.9278	0.8091	0.6930	0.8242
	sub-pre with WIIM (M4)	25.19	0.0081	5,853	0.9429	0.8510	0.7123	0.8821
6	whole-pre without WIIM (M1)	256.42	0.0028	4,267	0.9443	0.8901	0.6576	0.8951
	sub-pre without WIIM (M2)	178.29	0.0026	4,978	0.9443	0.9218	0.7005	0.9241
	whole-pre with WIIM (M3)	42.76	0.00089	7,774	0.9528	0.9091	0.7212	0.9258
	sub-pre with WIIM (M4)	25.19	0.00089	4,923	0.9619	0.9315	0.7406	0.9367

Method	Condition number	Relative AUC	Residual norm	Number of Iteration	nR	nD	nSAD	nSSD
non-preconditioning	648,496	100%	1.1479	3,963	0.5440	0.6517	0.4562	0.4809
whole-pre without WIIM	268.91	34.79%	1.1426	7,910	0.9079	0.7762	0.6106	0.7902
sub-pre without WIIM	189.75	45.27%	1.1437	5,792	0.9270	0.8280	0.6513	0.8593
whole-pre with WIIM	49.08	34.79%	0.0081	8,085	0.9173	0.8075	0.6898	0.8217
sub-pre with WIIM	26.85	45.27%	0.0081	5,853	0.9379	0.8478	0.7065	0.8712

Improving mesoscopic fluorescence molecular tomography via preconditioning and regularization

Abstract

Corrections

1. Introduction

2. Method

2.1 Configuration of the MFMT optical system

2.2 Forward model and inverse problem

2.3 Preconditioning to reduce the coherence of the sensitivity matrix

2.4 WIIM to speed-up the convergence rate of the inverse solver

3. Evaluation metrics for an objective assessment of the proposed methods

3.1 The effect of the condition number on the precision of the solution

3.2 Coherence of the sensitivity matrix

3.3 Signal to Noise ratio

3.4 Quantitative metrics to assess 3D reconstructions

4. In silico and experimental validation

4.1 Simulation settings

4.2 Coherence reduction via preconditioning

4.3 Faster convergence via the weighted iteration improvement method

4.4 Performance of preconditioning methods on experimental data

5. Conclusion

Funding

Disclosures

References and links

Cited By

Figures (7)

Tables (2)

Equations (22)

Biomedical Optics Express