Nonlinear greedy sparsity-constrained algorithm for direct reconstruction of fluorescence molecular lifetime tomography

Chuangjian Cai; Lin Zhang; Wenjuan Cai; Dong Zhang; Yanlu Lv; Jianwen Luo

doi:10.1364/BOE.7.001210

1. Introduction

Fluorescence molecular tomography (FMT) is increasingly utilized for gene therapy, cancer diagnosis and drug discovery [1, 2 ]. According to the illumination source, FMT can be classified into three modes: the continuous wave (CW) mode, the frequency-domain (FD) mode, and the time-domain (TD) mode. Compared to the other modes, the TD mode can acquire the temporal point spread functions and provide information at all frequencies, so it can be used to reconstruct fluorescence molecular lifetime tomography (FMLT). Fluorescence molecular lifetime is sensitive to physiological factors of the environment, such as tissue oxygenation, pH and glucose concentration [3]. Fluorescence molecular lifetime has been adopted to observe the fluorescence response energy transfer [4], to reduce the crosstalk between different fluorescent targets [5] and to improve the resolution of molecular yield tomography [6].

As for the reconstruction of FMLT, there are mainly two categories of methods. One is based on transforms, such as Fourier transform [7] and Laplace transform [8], which utilize only part of the characteristic information and rely on the selection of transformation factors. The other is directly based on the TD strategy and can make the best of all the measurement data. We previously developed an L1 regularization projected steepest descent (L1PSD) algorithm based on the TD strategy [9]. Simulations and phantom experiments demonstrated that the L1PSD algorithm performed well in localization of fluorescent targets and quantification of fluorescence lifetime, and provided high contrast between different fluorescent targets. However, the main drawbacks of L1PSD are that its resolution is poor and its solutions are not sparse enough for small targets. In the reconstructed FMLT images, the true targets are submerged in strong background noises. For the cancerous tumors in the early stages, the targets labeled with the fluorescent probes are typically small, namely sparse [10]. So effective sparsity-constrained methods for FMLT reconstruction are of importance. In order to enhance the spatial resolution of FMLT and obtain sparse reconstruction results for small-target problems, a method that gives higher priority to the values of the true targets (i.e., effective values) is needed. Greedy methods are known to perform well in this respect [11–13 ]. They can find an approximate solution of L₀ minimization problem, and are more suitable for problems with very sparse solution [14].

Compressed sensing (CS) is effective for sparsity-constrained reconstruction, both theoretically and practically [15–21 ]. Based on the CS framework, even far fewer measurements than those required by the Shannon sampling theory are needed [22, 23 ]. In the field of medical imaging, CS has been successfully used in magnetic resonance imaging [22], computed tomography [24], ultrasound imaging [25, 26 ] and photoacoustic tomography [27]. However, most CS studies are based on the assumption of a linear measurement model, which is not suitable for FMLT. The reconstruction of FMLT is a high-dimensional and ill-posed nonlinear problem. When the targets are sparse, nonlinear compressed sensing (NCS), of which the measurement of samples is nonlinear, can be employed in the reconstruction. For NCS, there are mainly two kinds of recovering methods, i.e., methods with L1 regularization [28–31 ] and greedy methods [32–36 ]. NCS has caught researchers’ attention these years [28–38 ]. However, which kind of recovering methods is better has not reached a consensus.

For NCS, L1 regularization is often adopted to induce sparsity in the solution, based on the assumption that the cost function is convex everywhere. L1 regularization is effective for regression with nonlinear models [28–30 ]. However, L1 regularization methods have some intrinsic drawbacks, which might cause unsatisfactory sparsity of the solution. Dobson et al. applied L1 regularization methods to the generalized linear models (GLM) [37]. In order to guarantee the error bound, namely the distance between the estimated and true statistical parameters, they introduced the notion of Restricted Strong Convexity (RSC). Around the true parameter, RSC requires a lower bound on the curvature of the cost function in a restricted set of directions. With a projected gradient descent method to solve the L1 regularized cost function, Agarwal et al. [31] guaranteed the accuracy by using a slightly different notion of RSC and restricted smoothness (RSM). Around the true optimum solution, RSM bounds the curvature of the cost function in a restricted set of directions. However, curvature bounds depend on the location of the true optimum solution for some cost functions, like loss functions in GLM [36]. Moreover, for the majority of L1 regularization methods, the true parameter is assumed to be bounded. For example, Bunea et al. recognized that the error bound for logistic regression with L1 regularization relies on the true parameter [28]. The selection of the regularization coefficient, based on the results reported in [29], implicitly requires the length of the true parameter to be short enough. So the methods with L1 regularization might not be suitable for all nonlinear models.

Greedy methods give higher priority to the most effective values and can find an approximate solution of L₀ minimization problem. They can perform well for nonlinear problems. Given some properties of the objective functions, such as the restricted diagonal deviation property, orthogonal matching pursuit (OMP) can find the exact support set for nonlinear models [32, 33 ]. Iterative hard thresholding algorithm can accurately recover sparse signals from few nonlinear observations, under conditions similar to linear compressed sensing [34]. Gradient hard thresholding pursuit algorithm performs well in sparse logistic regression and sparse precision matrix estimation [35]. Gradient support pursuit algorithm can guarantee a bounded distance between the true sparse optimum solution and the reconstructed sparse vector when the cost function has a stable restricted hessian matrix or a stable restricted linearization [36]. Therefore, greedy methods can be powerful for sparsity-constrained problems.

In order to improve the resolution of FMLT reconstruction, greedy methods are adopted in this paper. Among different greedy methods, OMP is widespread for its simplicity and high efficiency. Nonlinear OMP is recommended for its precise reconstruction given some properties of the objective functions, such as restricted diagonal deviation property [32, 33 ]. In this paper, an algorithm based on nonlinear OMP is proposed for FMLT reconstruction. The proposed algorithm is named accelerated nonlinear orthogonal matching pursuit (ANOMP) algorithm. It is divided into two parts, i.e., the outer iterations and the inner iterations. Each outer iteration selects multiple elements, instead of only one, to expand the support set of the inverse lifetime $I_{τ}$ , so as to accelerate the convergence. The inner iterations obtain an intermediate estimate of $I_{τ}$ , based on the support set estimated in the outer iterations. The selection of elements in the outer iterations is based on the gradients of mismatch error. For inner iterations, L1PSD is adopted to guarantee the stability and acquire accurate estimation results [9]. The stopping criterion for the outer iterations is based on the stability of the maximum reconstructed inverse lifetime values and is robust for problems with targets at different edge-to-edge distances (EEDs). ANOMP is based on the TD strategy directly, so it can make the best of all the data collected. Phantom experiments under different conditions and in vivo mouse experiments are performed to demonstrate that ANOMP can provide high resolution, high accuracy of localization and quantification, and good contrast at the same time.

2. Methods

2.1 Forward model for light propagation

To calculate the distribution of the excitation light and the Green’s function of the emission light, the telegraph equation [39] is selected as the forward model,

\begin{array}{l} \frac{D (r)}{c^{2}} \frac{\partial Φ^{2} (r, t)}{\partial t^{2}} + \frac{1}{c} (3 D (r) μ_{a} + 1) \frac{\partial Φ (r, t)}{\partial t} + μ_{a} Φ (r, t) \\ - \nabla \cdot [D (r) \nabla Φ (r, t)] = S (r, t), \end{array}

where

D (r)

denotes the diffusion coefficient, c denotes the velocity of light, r is the spatial coordinate, t is the time,

Φ (r, t)

denotes the density of photon and

S (r, t)

is the source term.

D (r) = 1 / (3 (μ_{a} + μ_{s}'))

, where

μ_{a}

and

μ_{s}'

are the absorption and reduced scattering coefficients. The standard Galerkin-finite element method (Galerkin-FEM) [39, 40 ] is adopted to transform Eq. (1) into a temporal recursion matrix equation. Then an implicit finite difference scheme [39, 41 ] is used for the approximation of time derivatives.

To reduce the difficulty of modeling, point source is used as the light source. $S_{x} (r, r_{s}, t) = δ (r - r_{s}, t)$ is the source term for the excitation light. For the emission light, the source term is determined by the distribution of the excitation light and fluorophore $E (r, t) = c η μ_{a f} (r) \exp (- t / τ (r)) / τ (r) = c η μ_{a f} (r) I_{τ} (r) \times$ $\exp (- I_{τ} (r) t)$ [8], where $η μ_{a f} (r)$ , $τ (r)$ and $I_{τ} (r)$ are the fluorescence yield, fluorescence lifetime and the inverse lifetime, respectively. So $S_{m} (r, t) = Φ_{x} (r, t) * E (r, t)$ [8], where * denotes the temporal convolution operator. The introduction of the inverse lifetime $I_{τ}$ is to prevent the singularity of the zero points in lifetime images [9]. For the background, $I_{τ}$ is assumed to be 0.

The emission signal detected can be formulated as

Φ_{m} (r_{s}, r_{d}, t) = \int_{Ω} Φ_{x} (r_{s}, r, t) * G_{m} (r, r_{d}, t) * E (r, t) d^{3} r,

where

G_{m} (r, r_{d}, t)

denotes the Green’s function of the emission light,

r_{d}

denotes the spatial coordinate of the detector and

r_{s}

denotes the spatial coordinate of the source.

2.2 FMLT reconstruction based on ANOMP algorithm

The image domain is discretized and decomposed into $N_{G}$ voxels, so the lengths of the yield map $η μ_{a f}$ , lifetime map $τ$ , and the inverse lifetime map $I_{τ}$ are all $N_{G}$ . With the pre-reconstructed yield map [9], the reconstruction of FMLT can be divided into the outer iterations and the inner iterations. The outer iterations expand the support set of $I_{τ}$ continuously. The elements with the highest gradients of L2 norm of the fitting error are considered to be effective and are selected. The inner iterations, i.e., the estimation steps of each outer iteration, obtain an intermediate estimate of $I_{τ}$ , based on the estimated support set.

To avoid dealing with the singularity of the zero points, the inverse lifetime map $I_{τ}$ , instead of FMLT image $τ$ , is reconstructed first [9]. With the yield map $η μ_{a f}$ obtained, the FMLT problem can be developed as

Φ_{m} = F (I_{τ}),

where the fluorescent measurement vector

Φ_{m}

is rearranged in order of excitation sources, detectors and time points, and has a length of

S \times D \times I

, where S is the number of excitation sources, D is the number of detectors and I is the number of discrete time points. The ith element of

Φ_{m}

and the nonlinear forward function F is given by

\begin{array}{l} Φ_{m}^{i} = F^{i} = \sum_{n = 1}^{N_{G}} F_{n}^{i} Δ V = \sum_{n = 1}^{N_{G}} Δ V [Φ_{x} (r_{s}^{i}, n, t) * G_{m} (n, r_{d}^{i}, t) * . \\ (c η μ_{a f} (n) I_{τ} (n) e^{- I_{τ} (n) t})] |_{t = t^{i}} \end{array}

Nonlinear OMP, as the algorithm of the outer iterations, is used to obtain a sparse solution of the FMLT problem (3).

OMP is an iterative greedy algorithm that finds an approximate solution of the L₀ minimization problem. Besides, for the target algorithm, L₁ regularization is adopted to obtain an intermediate estimate in the inner iterations. Similar to that of linear OMP, the optimization function of the nonlinear OMP for inverse lifetime is given by

\min \begin{matrix}  \end{matrix} | | I_{τ} | |_{0} s . t . \begin{matrix}  \end{matrix} F (I_{τ}) = Φ_{m} .

In each iteration step, new elements are selected,

i_{k} = \arg_{i} \max | \nabla F^{T} (I_{τ}^{k - 1}) (Φ_{m} - F (I_{τ}^{k - 1})) |,

where

\nabla F

is the local gradient of F. The selection criterion of elements is based on the gradient of mismatch error

0 .5|| (Φ_{m} - F (I_{τ}^{k - 1})) {||}_{2}^{2}

, which is referred to as effectiveness gradient in this paper.

Because FMLT is a high-dimensional nonlinear problem, multiple elements, instead of only one, are selected in each iteration to accelerate the convergence. So the selection criterion turns to be

z = | \nabla F^{T} (I_{τ}^{k - 1}) (Φ_{m} - F (I_{τ}^{k - 1})) |,

Z = S u p p (z_{N_{s}}),

where z is the magnitude of the effectiveness gradient. Equation (8) extracts the indices of

N_{S}

components of vector z that have the largest magnitudes. The gradient of F and the corresponding elements are formulated as

\nabla F = (\begin{matrix} \nabla F_{1}^{1} & \dots & \nabla F_{N_{G}}^{1} \\ ⋮ & ⋱ & ⋮ \\ \nabla F_{1}^{S \times D \times I} & \dots & \nabla F_{N_{G}}^{S \times D \times I} \end{matrix}),

\begin{array}{l} \nabla F_{n}^{i} = Δ V [Φ_{x} (r_{s}^{i}, n, t) * G_{m} (n, r_{d}^{i}, t) * . \\ (c η μ_{a f} (n) (1 - I_{τ} (n) t) e^{- I_{τ} (n) t})] |_{t = t^{i}} \end{array}

Then the new indices should be merged with the estimated support set, as

Γ_{k} = Γ_{k - 1} \cup Z,

where

Γ

is the support set.

With the updated support set, the inner iterations are conducted to obtain an intermediate estimate of $I_{τ}$ . According to the principle of traditional greedy methods [33], the new estimate should be $I_{τ}^{k} = \arg$ $\min 0.5 | | F (I_{τ}) - Φ_{m} | |_{2}^{2}$ s.t. $I_{τ} |_{Γ_{k}} \geq 0; I_{τ} |_{Γ_{k}^{c}} = 0$ , where c is the complement operator. However, the reconstruction of FMLT is a severely ill-posed problem, because of the huge amount of noise, the severe redundancy of measurement data, the strong scattering of light and the high nonlinearity of the forward model (Eq. (3)) under physical conditions. For reconstruction of fluorescence yield tomography, regularization methods such as the truncated singular value algorithm are introduced in greedy algorithms to reduce the ill-posedness [11]. However, the direct reconstruction of FMLT is a nonlinear problem, so the truncated singular value algorithm is not suitable. In this paper, L1PSD method is employed, which can guarantee the stability of convergence and enhance the anti-noise capability of sparsity-constrained reconstruction [9]. So the optimization function for inner iterations is given by

\min 0.5 (| | F (I_{τ}) - Φ_{m} | |_{2}^{2} + λ \sqrt{| | I_{τ} | |_{2}^{2} + μ}) s . t . I_{τ} |_{Γ_{k}} \geq 0; I_{τ} |_{Γ_{k}^{c}} = 0.

In order to reduce the computational complexity, the components of F in Eq. (3) and columns of $\nabla F$ corresponding to $Γ_{k}^{c}$ are removed. The number of the indices of $Γ_{k}$ is written as $N_{G}^{k}$ . Then $F^{i}$ and $\nabla F$ turn to be

\begin{array}{l} Φ_{m}^{i} = F^{i} = \sum_{n = 1}^{N_{G}^{k}} F_{Γ_{k}^{n}}^{i} Δ V = \sum_{n = 1}^{N_{G}^{k}} Δ V [Φ_{x} (r_{s}^{i}, Γ_{k}^{n}, t) * G_{m} (Γ_{k}^{n}, r_{d}^{i}, t) *, \\ (c η μ_{a f} (Γ_{k}^{n}) I_{τ} (Γ_{k}^{n}) e^{- I_{τ} (Γ_{k}^{n}) t})] |_{t = t^{i}} \end{array}

\nabla F = (\begin{matrix} \nabla F_{Γ_{k}^{1}}^{1} & \dots & \nabla F_{Γ_{k}^{N_{G}^{k}}}^{1} \\ ⋮ & ⋱ & ⋮ \\ \nabla F_{Γ_{k}^{1}}^{S \times D \times I} & \dots & \nabla F_{Γ_{k}^{N_{G}^{k}}}^{S \times D \times I} \end{matrix}),

where

Γ_{k}^{c}

denotes the ith element of the support set. Mathematically, there is no order for the support set. The elements of the support set are assumed to be arranged in ascending order in this study.

For the inner iterations, the stopping criterion is modified as $Δ R / R_{0} < ξ$ , where $R_{0}$ is the initial residual error, $R_{0} = \sqrt{| | F (i n n_I_{τ}^{0}) - Φ_{m} | |_{2}^{2}}$ , and $Δ R$ is the difference of residual errors of adjacent iterations, $Δ R = \sqrt{| | F (i n n_I_{τ}^{l}) - Φ_{m} | |_{2}^{2}} - \sqrt{| | F (i n n_I_{τ}^{l - 1}) - Φ_{m} | |_{2}^{2}}$ . The stopping criterion takes two kinds of situations into consideration. The first one is that the iterative process reaches a convergence for the inner iterations, for which $0 \leq Δ R / R_{0} < ξ$ . The second one is that there is an oscillation of convergence, for which $Δ R < 0$ . The second situation usually occurs at the very beginning of the outer iterations, when the support set is too small. When either of the two situations occurs, the stopping criterion takes effect and the inner iterations are interrupted.

For outer iterations, the elements possessing the highest effectiveness gradients are considered to be the most effective and are selected. The stopping criterion is the key factor that influences the sparsity and accuracy of the reconstruction results. If the true sparsity is known, the stopping criterion can be set as the lower limit of reconstruction sparsity. However, the true sparsity is difficult to obtain in practice. The stopping criterion can also be based on the effectiveness gradient, like the upper limit of the maximum effectiveness gradient, which means that the outer iterations are interrupted if the maximum effectiveness gradient exceeds a preset upper limit. However, the effectiveness gradients of the background are strongly influenced by the EED of targets, so the stopping criterion based on the effectiveness gradients is not robust.

For non-greedy methods, the stopping criteria can be based on the difference of reconstruction results of adjacent iterations [42]. In [8], the maximum reconstructed inverse lifetime values of the targets are confirmed to be good approximations of real values of targets and have great significance for quantification. For ANOMP, the difference of reconstruction results of different outer iterations is described by the difference of the maximum reconstructed inverse lifetime values. Figure 1(a) is the maximum reconstructed inverse lifetime values of different outer iterations of Expt. 1 described below, and Figs. 1(b)-1(e) are the corresponding reconstructed inverse lifetime maps of the 1st, 4th, 8th and 13th outer iterations, respectively. At the very beginning of the outer iterations, when the reconstruction sparsity is much smaller than the real one, the maximum reconstructed values are not close to the real values and are unstable, because the support set has not covered the real distribution of the targets and the reconstruction algorithm of the inner iterations is difficult to recover the real values. So the reconstruction with a too small support set cannot provide accurate and stable results. After the support set covers the real distribution of the targets, the reconstructed inverse lifetime values of the targets will remain close to the real values, so the maximum reconstructed values can keep stable. When the maximum reconstructed values remain stable, the iterative process is thought to reach a convergence and excessive iterations will only affect the resolution and sparsity, so the iterations are interrupted.

Fig. 1 (a) The maximum reconstructed inverse lifetime values of different outer iterations in Expt. 1. Horizontal axis denotes the number of outer iteration and vertical axis denotes the maximum inverse lifetime (unit: s⁻¹). (b-e) The reconstructed inverse lifetime maps of the 1st, 4th, 8th and 13th outer iterations (unit: s⁻¹), the maximum values of which are $0.11 \times 10^{9} s^{- 1}$ , $1.28 \times 10^{9} s^{- 1}$ , $1.08 \times 10^{9} s^{- 1}$ and $1.09 \times 10^{9} s^{- 1}$ , respectively.

Download Full Size | PDF

Set $Γ_{v}$ to be the set of maximum reconstructed values of the outer iterations that belong to the stable part. The range of $Γ_{v}$ should be below a preset lower limit, i.e.,

\max (Γ_{_{v}}^{k}) - \min (Γ_{_{v}}^{k}) < ω .

If the range of $Γ_{v}$ violates (15), $Γ_{v}$ will be reset, and only the maximum reconstructed value of the current step will be included, i.e., $Γ_{_{v}}^{k} = {\max (I_{τ}^{k})}$ .

The number of outer iterations is a compromise between accuracy and sparsity. In Fig. 1(a), the maximum reconstructed values of the 4th and 5th outer iteration are very close. However, the iterations only reach a pseudo optimum solution and the maximum reconstructed values are still unstable. In order to avoid the pseudo optimum solution, the maximum reconstructed values of enough number of adjacent outer iterations should be close to each other. So a stopping criterion is formulated as

\frac{\max (Γ_{_{v}}^{k}) - \min (Γ_{_{v}}^{k}) + ζ}{n u m (Γ_{_{v}}^{k}) - 1} < ζ,

where the operator

n u m

obtains the number of elements of

Γ_{_{v}}^{k}

. When the range of

Γ_{_{v}}^{k}

,

\max (Γ_{_{v}}^{k}) - \min (Γ_{_{v}}^{k})

, is small enough and the number of the outer iterations that belong to the stable part,

n u m (Γ_{_{v}}^{k})

, is large enough, the outer iterations are interrupted. The term

ζ

added in the numerator part of (16) is to deal with the situation when only two elements constitute

Γ_{_{v}}^{k}

and their difference is extremely small. When the range of

Γ_{_{v}}^{k}

is relatively small, the number of elements of

Γ_{_{v}}^{k}

and the number of outer iterations are allowed to be small, which can lead to good reconstruction sparsity.

The main steps of the proposed algorithm are summarized as below.

Algorithm 1 Accelerated Nonlinear Orthogonal Matching Pursuit (ANOMP)

Parameters setting

For the outer iterations, set the number of the selected elements of each outer iteration , the lower limit of stopping criterion $ζ$ , the lower limit of the range of the maximum reconstructed values $ω$ , the initial inverse lifetime $I_{τ}^{0} = 0$ , the initial iteration number k = 1, the initial support set $Γ_{0} = Ø$ ; for the inner iterations, set the initial inverse lifetime $i n n I_{τ}^{0} = 10^{8} s^{- 1}$ , the initial regularization parameter $λ_{0}$ , the attenuation coefficient of the regularization parameter $D_{e} =0 .6$ , the maximum ratio between the difference of residual errors of adjacent iterations and the initial residual error $ξ$ .

ANOMP

(1) Compute the local gradient $z = | \nabla F^{T} (I_{τ}^{k - 1}) (Φ_{m} - F (I_{τ}^{k - 1})) |$ ;
(2) Extract the indices of the elements that are the most effective: $z = | \nabla F^{T} (I_{τ}^{k - 1}) (Φ_{m} - F (I_{τ}^{k - 1})) |$ , $Z = S u p p (z_{N_{s}})$ ;
(3) Update the support set $Γ_{k} = Γ_{k - 1} \cup Z$ ;
(4) Conduct the inner iterations over the updated support set.
(4.1) The initial inverse lifetime $i n n I_{τ}^{0} = 10^{8} s^{- 1}$ , the length of which is equal to the number of elements of $Γ_{k}$ , namely $N_{G}^{k}$ . The components of F and columns of $\nabla F$ corresponding to $Γ_{k}^{c}$ are removed. The original search direction $d_{0} = - {(\nabla F (i n n I_{τ}^{0}))}^{T} (F (i n n I_{τ}^{0}) - Φ_{m}) + λ i n n I_{τ}^{0} / (2 \sqrt{| | i n n I_{τ}^{0} | |_{2}^{2} + μ})$ , the original step size $t_{0} = | | d_{0} | |_{2}^{2} / (| | \nabla F (i n n I_{τ}^{0}) d_{0} | |_{2}^{2} + λ_{0} | | d_{0} | |_{2}^{2} / (2 \sqrt{|| i n n I_{τ}^{0} {||}_{2}^{2} + μ}))$ and the initial residual error is $R_{0} = \sqrt{| | F (i n n I_{τ}^{0}) - Φ_{m} | |_{2}^{2}}$ . $l = l + 1$ ;
(4.2) Update $i n n I_{τ}^{l}$ with step size $t_{l - 1}$ and search direction $d_{l - 1}$ , $i n n I_{τ}^{l} = i n n I_{τ}^{l - 1} + t_{l - 1} d_{l - 1}$ , and set the negative part of $i n n I_{τ}^{l}$ to be zero. Update $λ_{l} = D_{e} λ_{l - 1}$ and $Δ R = \sqrt{| | F (i n n I_{τ}^{l}) - Φ_{m} | |_{2}^{2}} -$ $\sqrt{| | F (i n n I_{τ}^{l - 1}) - Φ_{m} | |_{2}^{2}}$ ;
(4.3) $l = l + 1$ . If $Δ R / R_{0} \geq ξ$ go to step 4.2;
(4.4) Obtain $I_{τ}^{k}$ corresponding to the reconstructed $i n n I_{τ}$ ;
(5) If k = 1, $Γ_{_{v}}^{k} = {\max {I_{τ}^{k}}}$ , go to step 7;
(6) $Γ_{_{v}}^{k} = Γ_{_{v}}^{k - 1} \cup {\max {I_{τ}^{k}}}$ . If $\max {Γ_{_{v}}^{k}} - \min {Γ_{_{v}}^{k}} \geq ω$ , $Γ_{_{v}}^{k} = {\max {I_{τ}^{k}}}$ ; else if $(\max {Γ_{_{v}}^{k}} - \min {Γ_{_{v}}^{k}} + ζ) / n u m {Γ_{_{v}}^{k}} < ζ$ , go to step 8;
(7) k = k + 1. Go to step 1;
(8) Obtain FMLT corresponding to the reconstructed $I_{τ}$ .

To compromise between the computational time and the robustness for convergence, the number of the selected elements of each outer iteration $N_{S}$ is empirically set to be 10. The lower limit of the stopping criterion $ζ$ and the lower limit of the range of the maximum reconstructed values $ω$ are set to be $0.01 \times 10^{9} s^{- 1}$ and $0.1 \times 10^{9} s^{- 1}$ for all the experiments, to ensure convergence of the outer iterations. For the inner iterations, the maximum ratio between the difference of residual errors of adjacent iterations and the initial residual error $ξ$ is set to be 0.007. After being normalized by the sum of the measured fluorescent values, the empirical range of the initial regularization parameter $λ_{0}$ is between $5 \times 10^{- 15}$ and $5 \times 10^{- 14}$ . If $λ_{0}$ is beyond $5 \times 10^{- 14}$ , the residual errors in the inner iterations will oscillate, no matter how many elements are included in the estimated support set. If $λ_{0}$ is below $5 \times 10^{- 15}$ , the reconstruction problem will turn into a severely ill-posed one and the reconstruction accuracy will decrease. In this paper, $λ_{0}$ is set to be $1 \times 10^{- 14}$ for all the experiments.

2.3 Quantitative metrics for FMLT reconstruction

In order to evaluate the performance of ANOMP, several quantitative metrics are selected.

The mean absolute error (MAE) measures the closeness between the reconstructed results and the real values. Because a small inverse lifetime value will turn into a very big lifetime value and the theoretical lifetime of the background is infinite, the MAE is suitable for the reconstructed inverse lifetime tomographic images, instead of the final lifetime images. The MAE is defined as

M A E = \frac{| | I_{τ}^{r e c o n} - I_{τ}^{t r u e} | |_{1}}{N_{G}},

where

I_{τ}^{r e c o n}

and

I_{τ}^{t r u e}

are the reconstructed inverse lifetime and the real inverse lifetime of the whole object, including the background and the targets.

The contrast-to-noise ratio (CNR) [43] is adopted to indicate if the real targets are well recovered or submerged in the noise, i.e., target detectability. The CNR is also suitable for the inverse lifetime and is defined as

C N R = \frac{\bar{I_{τ}^{R O I}} - \bar{I_{τ}^{B C K}}}{{(w_{R O I} σ_{R O I}^{2} + w_{B C K} σ_{B C K}^{2})}^{1 / 2}},

where

\bar{I_{τ}^{R O I}}

and

\bar{I_{τ}^{B C K}}

are the means of the reconstructed inverse lifetime values in the region-of-interested (ROI), i.e., the region of the real targets, and the background,,

w_{R O I}

and

w_{B C K}

are the volumes of the ROI and background (relative to the whole phantom),

σ_{R O I}^{2}

and

σ_{B C K}^{2}

are the variances of the inverse lifetime in the ROI and background, respectively.

To evaluate the quantification accuracy of lifetime, the relative error (RE) is defined as

R E = \frac{| τ^{r e c o n} - τ^{t r u e} |}{τ^{t r u e}},

where

τ^{r e c o n}

and

τ^{t r u e}

are the reconstructed and real lifetime values of a fluorescent target, respectively.

τ^{r e c o n}

is defined as the valley lifetime value of the reconstructed target. In this paper, the bigger value of the REs of two fluorescent targets is selected to be an evaluation metric, called as the maximum relative error (

R E_{\max}

).

In order to evaluate the resolution of the reconstructed lifetime images, the separability of two fluorescent targets is considered. If the lifetime values along the line between the centers of the two targets are all smaller than 3 times of the minimum reconstructed lifetime value of the two targets, the targets are considered to be inseparable; or otherwise the targets are separable.

All the evaluation metrics are used in the phantom experiments. In the in vivo animal experiments, the exact locations of targets are unknown, so only $R E_{\max}$ and the separability are evaluated.

3. Experiments and results

3.1 Phantom experiments

In order to test the performance of ANOMP, three phantom experiments, Expt. 1, Expt. 2 and Expt. 3, with fluorophores at different EEDs (9 mm, 5 mm and 2 mm), were conducted. A fiber-coupled, time-correlated and single-photon-counting (TCSPC) system was built, the schematic of which could be found in our previous work [9]. A femtosecond laser generator (Spectra-Physics, Newport Corporation, Irvine, CA) provided the repetitive ultrashort 775 nm-wavelength pulse (80 MHz, 100 fs pulse-width). The detector consisted of a multi–anode photomultiplier (PMT) (PML-16-C, Becker and Hickl, Berlin, Germany) and the TCSPC module. Different from the system previously developed [9]), a group of filters, including an achromatic doublet (AC254-030-B, Thorlabs, Newton, NJ), two bandpass fluorescence filters with a center wavelength of 840 nm (ff01-840/12-25, Semrock, Rochester, NY; XBPA840, Asahi Spectra, Torrance, CA) and a long-pass edge filter (BLP01-785R-12.5, Semrock, Rochester, NY), was set in front of the detector, to eliminate most of the excitation light . The measurement data were subsampled and truncated, to reduce the computational cost. The final temporal resolution was set to be 75 ps and time gate was set to be 5.175 ns. A set of measurement data consisted of 18 projections, and each projection included 11 probe points and a field of view (FOV) of 220°. So S = 18, D = 11, and I = 69. As for the phantom setup, a cylindrical phantom (30 mm in diameter) filled with 1% intralipid solution (with absorption coefficient of 0.02 cm⁻¹ and reduced scattering coefficient of 10 cm⁻¹) was used as the background. Two small cylinders (4 mm in diameter and 10 mm in height) filled with 10 μM indocyanine green/dimethyl sulphoxide (ICG/DMSO) were used as fluorescent targets. The real lifetime and inverse lifetime of the targets are $0 .97 \times 10^{- 9} s$ and $1 .03 \times 10^{9} s^{- 1}$ [44]. The ICG has a peak excitation wavelength at 780 nm in DMSO and a peak emission wavelength at 830 nm, which are close to the excitation wavelength and the central wavelength of the fluorescence filters in the system [45].

In the reconstruction, only the part containing fluorophores was taken into account. The target domain was 30 mm in diameter and 10 mm in height. The geometry model was discretized into 3,145 nodes and 15,407 tetrahedral elements. The Galerkin-FEM [39, 40 ] was adopted to calculate the distribution of the excitation light and the Green’s function of the emitted light. Then the target domain was discretized into 4,182 voxels of size $1 \times 1 \times 2 m m^{3}$ , and the distribution of excitation light and Green’s function of the emitted light were also interpolated.

Figures 2, 3 and 4 are the reconstruction results of Expt. 1, Expt. 2 and Expt. 3, respectively. The reconstruction results of L1PSD [9] are also presented for comparison. Figures 2(a), 2(b), 3(a), 3(b), 4(a) and 4(b) are the reconstructed inverse lifetime tomographic images of ANOMP and L1PSD, respectively. As can be seen, the reconstructed inverse lifetime values of the background using L1PSD (Figs. 2(b), 3(b) and 4(b) ) are too big. Especially in Figs. 4(b), with an EED of 2 mm, the reconstructed inverse lifetime values of the background are almost the same as those of the targets. The resolutions of ANOMP (Figs. 2(a), 3(a) and 4(a) ) are much higher, and the two targets can be separated thoroughly. Figures 2(c), 3(c) and 4(c) show the inverse lifetime profiles along the dotted lines indicated in the reconstructed inverse lifetime tomographic images. The red lines, the blue lines and the green lines correspond to the real values, the results of ANOMP and L1PSD, respectively. When the EEDs of two targets are not so small (Figs. 2(c) and 3(c) ), the peak reconstructed values of ANOMP and L1PSD are close to the real values. However, when two targets are close (Fig. 4(c)), the strong reconstructed background of L1PSD can even affect the quantification accuracy. Figures 2(d), 2(e), 3(d), 3(e), 4(d) and 4(e) are the corresponding FMLT images. Because the reconstructed inverse lifetime values close to 0 will turn into extremely big lifetime values, the lifetime values exceeding 1/20% = 5 times of the real lifetime value are eliminated in the FMLT images. The minimum values of the color bars for the FMLT images in Figs. 2, 3 and 4 are set to be the corresponding minimum value of the images, respectively. Figures 2(f), 3(f) and 4(f) show the lifetime profiles along the dotted lines indicated in the FMLT images.

Fig. 2 The reconstruction results of Expt. 1 (EED = 9 mm). (a) and (b) are the reconstructed inverse lifetime tomographic images of ANOMP and L1PSD (unit: s⁻¹). (c) shows the inverse lifetime profiles along the dotted lines indicated in (a) and (b). Horizontal and vertical axes denote location (unit: mm) and inverse lifetime (unit: s⁻¹), respectively. (d) and (e) are the FMLT images (unit: s) corresponding to (a) and (b). (f) shows the lifetime profiles along the dotted lines indicated in (d) and (e). Horizontal and vertical axes denote location (unit: mm) and lifetime (unit: s), respectively. In (c) and (f), the red, blue and green lines correspond to the real values, the results of ANOMP and L1PSD, respectively.

Download Full Size | PDF

Fig. 3 The reconstruction results of Expt. 2 (EED = 5 mm). (a) and (b) are the reconstructed inverse lifetime tomographic images of ANOMP and L1PSD (unit: s⁻¹). (c) shows the inverse lifetime profiles along the dotted lines indicated in (a) and (b). Horizontal and vertical axes denote location (unit: mm) and inverse lifetime (unit: s⁻¹), respectively. (d) and (e) are the FMLT images (unit: s) corresponding to (a) and (b). (f) shows the lifetime profiles along the dotted lines indicated in (d) and (e). Horizontal and vertical axes denote location (unit: mm) and lifetime (unit: s), respectively. In (c) and (f), the red, blue and green lines correspond to the real values, the results of ANOMP and L1PSD, respectively.

Download Full Size | PDF

Fig. 4 The reconstruction results of Expt. 3 (EED = 2 mm). (a) and (b) are the reconstructed inverse lifetime tomographic images of ANOMP and L1PSD (unit: s⁻¹). (c) shows the inverse lifetime profiles along the dotted lines indicated in (a) and (b). Horizontal and vertical axes denote location (unit: mm) and inverse lifetime (unit: s⁻¹), respectively. (d) and (e) are the FMLT images (unit: s) corresponding to (a) and (b). (f) shows the lifetime profiles along the dotted lines indicated in (d) and (e). Horizontal and vertical axes denote location (unit: mm) and lifetime (unit: s), respectively. In (c) and (f), the red, blue and green lines correspond to the real values, the results of ANOMP and L1PSD, respectively.

Download Full Size | PDF

Table 1 lists the quantitative metrics of the reconstruction results of the phantom experiments. The MAEs of ANOMP are all lower than 1/3 of those of L1PSD. The lower MAEs of ANOMP indicate that ANOMP has higher quantification accuracy than L1PSD. The higher CNRs of ANOMP demonstrate that the reconstructed targets can be detected better with ANOMP. For Expt. 1 and Expt. 2, L1PSD can guarantee the quantification accuracy of the reconstructed targets and has low $R E_{\max}$ . However, for Expt. 3, when the EED is only 2 mm and the reconstructed background is too strong, L1PSD has a large $R E_{\max}$ (27.48%). ANOMP can provide high quantification accuracy of the reconstructed targets with different EEDs ( $R E_{\max}$ <10%). As for the separability, the two targets cannot be separated with L1PSD. By contrast, ANOMP can effectively separate the two targets.

Table 1. Quantitative metrics of reconstruction results of the phantom experiments

View Table

3.2 In vivo mouse experiments

To further test the performance of ANOMP in practical applications, in vivo mouse experiments were conducted. The experiments were approved by the Ethics Committee of Tsinghua University. An 8-week BALB/c nude mouse was anesthetized by 0.2 ml 1% Nembutal during the experiments. Two small cylinders (4 mm in diameter and 10 mm in height) filled with 10 μM ICG/DMSO were used as the fluorescent targets, and were implanted into the abdomen of the mouse. The abdomen could be approximated with homogeneous model [46]. Then the mouse was put into a cylindrical container (30 mm in diameter) and its arms and legs were stuck to the inside wall of the container with waterproof rubberized fabric. 1% intralipid solution, with optical parameters close to the abdomen of mouse [47], was poured into the container and covered the abdomen of the mouse. The methods and parameters of data collection and reconstruction were the same as those used in the phantom experiments.

The reconstruction results are shown in Fig. 5 . Figures 5(a) and 5(b) are the reconstructed inverse lifetime tomographic images of ANOMP and L1PSD. Figure 5(c) shows the inverse lifetime profiles along the dotted lines indicated in Figs. 5(a) and (b). The reconstructed targets of ANOMP are separated thoroughly, while those of L1PSD are submerged in the strong background. Figures 5(d) and 5(e) are the corresponding FMLT images. Figure 5(f) shows the lifetime profiles along the dotted lines indicated in Figs. 5(d) and 5(e). The exact locations of the targets are difficult to acquire with the single-mode system, so the quantification evaluation of the reconstruction results only adopts the parameters of $R E_{\max}$ and separability in the in vivo experiments. The $R E_{\max}$ of ANOMP and L1PSD are as small as 3.65% and 4.10%, respectively. The reconstructed targets of ANOMP are separable, while those of L1PSD are inseparable.

Fig. 5 The reconstruction results of the in vivo experiment. (a) and (b) are the reconstructed inverse lifetime tomographic images of ANOMP and L1PSD (unit: s⁻¹). (c) shows the inverse lifetime profiles along the dotted lines indicated in (a) and (b). Horizontal and vertical axes denote location (unit: mm) and inverse lifetime (unit: s⁻¹), respectively. (d) and (e) are the FMLT images (unit: s) corresponding to (a) and (b). (f) shows the lifetime profiles along the dotted lines indicated in (d) and (e). Horizontal and vertical axes denote location (unit: mm) and lifetime (unit: s), respectively. In (c) and (f), the blue lines and green lines correspond to the reconstruction results of ANOMP and L1PSD, respectively.

Download Full Size | PDF

4. Discussions

For FMLT reconstruction, the L1PSD algorithm based on the TD strategy has previously been demonstrated to perform well in localization of fluorescent targets and quantification of fluorescence lifetime, and provide good contrast between different fluorescence targets [9]. However, L1PSD has poor resolution and its solution is not sparse enough for small targets. Moreover, when the EED of two targets is small (Expt. 3), the strong background noise could even affect the quantification accuracy of the reconstructed targets ( $R E_{\max}$ = 27.48%). In order to improve the resolution of FMLT, ANOMP, a nonlinear greedy sparsity-constrained algorithm, is proposed. Greedy methods can find an approximate solution of L₀ minimization problem. ANOMP consists of two parts, i.e., the outer iterations and the inner iterations. The outer iterations expand the support set continuously, based on the gradients of mismatch error. The inner iterations utilize the L1PSD algorithm to obtain an intermediate estimate of $I_{τ}$ , based on the estimated support set. The stopping criterion for the outer iterations is based on the stability of the maximum reconstructed values and is robust for problems with targets at different EEDs. Only after the support set covers the real distribution of the targets, the reconstructed inverse lifetime values of the targets can be accurate and remain close to the real values; and the maximum reconstructed values can keep stable.

Phantom experiments with two fluorophores at different EEDs were carried out to test the performance of ANOMP. The reconstruction results indicate that ANOMP can provide high quantification accuracy, good detectability of the fluorescent targets, and good separability of two targets. Especially when the EED is only 2 mm, the reconstructed background is too strong and thus affects the quantification accuracy for L1PSD. In comparison, ANOMP still behaves well in quantification of the fluorescent targets. In addition, in vivo mouse experiments were conducted to further evaluate the performance of ANOMP in practical applications. ANOMP still provides high quantification accuracy for the fluorescent targets and thoroughly separates the fluorophores implanted in the abdomen of the mouse. By contrast, the reconstructed background of L1PSD is too strong and makes the reconstructed targets inseparable. The results of the phantom experiments and in vivo mouse experiments demonstrate that, ANOMP is effective for the problems with sparse targets.

5. Conclusion

In conclusion, an accelerated nonlinear orthogonal matching pursuit algorithm, as a kind of nonlinear greedy sparsity-constrained reconstruction method, is proposed to enhance the spatial resolution and sparsity of the solution for fluorescence molecular lifetime tomography problems with sparse targets. The reconstruction results of the phantom experiments with different edge-to-edge distances and the in vivo mouse experiments demonstrate that ANOMP can provide high quantification accuracy for the fluorescent targets, even if the EED is small, and high resolution. However, the selection of regularization parameter, i.e., λ0, is still empirical. Future work will be focused on finding a formula for the determination of λ0, improving the shapes of the reconstructed targets and applying FMLT to in vivo cases.

Acknowledgments

This work is supported by the National Natural Science Foundation of China under Grant Nos. 81227901, 81271617, 61322101 and 61361160418, and the National Major Scientific Instrument and Equipment Development Project under Grant No. 2011YQ030114. We thank Profs. Hongen Liao and Sen Song in the Department of Biomedical Engineering at Tsinghua University for providing optical fibers and Nembutal, respectively. We also thank Huizhi Guang in our laboratory for helpful discussions.

References and links

1. V. Ntziachristos, J. Ripoll, L. V. Wang, and R. Weissleder, “Looking and listening to light: the evolution of whole-body photonic imaging,” Nat. Biotechnol. 23(3), 313–320 (2005). [CrossRef] [PubMed]

2. J. K. Willmann, N. van Bruggen, L. M. Dinkelborg, and S. S. Gambhir, “Molecular imaging in drug development,” Nat. Rev. Drug Discov. 7(7), 591–607 (2008). [CrossRef] [PubMed]

3. M. A. O’Leary, D. A. Boas, X. D. Li, B. Chance, and A. G. Yodh, “Fluorescence lifetime imaging in turbid media,” Opt. Lett. 21(2), 158–160 (1996). [CrossRef] [PubMed]

4. J. McGinty, D. W. Stuckey, V. Y. Soloviev, R. Laine, M. Wylezinska-Arridge, D. J. Wells, S. R. Arridge, P. M. French, J. V. Hajnal, and A. Sardini, “In vivo fluorescence lifetime tomography of a FRET probe expressed in mouse,” Biomed. Opt. Express 2(7), 1907–1917 (2011). [CrossRef] [PubMed]

5. S. S. Hou, W. L. Rice, B. J. Bacskai, and A. T. N. Kumar, “Tomographic lifetime imaging using combined early- and late-arriving photons,” Opt. Lett. 39(5), 1165–1168 (2014). [CrossRef] [PubMed]

6. W. L. Rice, S. Hou, and A. T. N. Kumar, “Resolution below the point spread function for diffuse optical imaging using fluorescence lifetime multiplexing,” Opt. Lett. 38(12), 2038–2040 (2013). [CrossRef] [PubMed]

7. R. E. Nothdurft, S. V. Patwardhan, W. Akers, Y. Ye, S. Achilefu, and J. P. Culver, “In vivo fluorescence lifetime tomography,” J. Biomed. Opt. 14(2), 024004 (2009). [CrossRef] [PubMed]

8. L. Zhang, F. Gao, H. He, and H. Zhao, “Three-dimensional scheme for time-domain fluorescence molecular tomography based on Laplace transforms with noise-robust factors,” Opt. Express 16(10), 7214–7223 (2008). [CrossRef] [PubMed]

9. C. Cai, L. Zhang, J. Zhang, J. Bai, and J. Luo, “Direct reconstruction method for time-domain fluorescence molecular lifetime tomography,” Opt. Lett. 40(17), 4038–4041 (2015). [CrossRef] [PubMed]

10. P. Mohajerani, A. A. Eftekhar, J. Huang, and A. Adibi, “Optimal sparse solution for fluorescent diffuse optical tomography: theory and phantom experimental results,” Appl. Opt. 46(10), 1679–1685 (2007). [CrossRef] [PubMed]

11. J. Shi, X. Cao, F. Liu, B. Zhang, J. Luo, and J. Bai, “Greedy reconstruction algorithm for fluorescence molecular tomography by means of truncated singular value decomposition conversion,” J. Opt. Soc. Am. A 30(3), 437–447 (2013). [CrossRef] [PubMed]

12. J. Yu, F. Liu, J. Wu, L. Jiao, and X. He, “Fast source reconstruction for bioluminescence tomography based on sparse regularization,” IEEE Trans. Biomed. Eng. 57(10), 2583–2586 (2010). [CrossRef] [PubMed]

13. 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]

14. J. Tropp and S. J. Wright, “Computational methods for sparse solution of linear inverse problems,” Proc. IEEE 98(6), 948–958 (2010). [CrossRef]

15. S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” SIAM J. Sci. Comput. 20(1), 33–61 (1998). [CrossRef]

16. M. A. Figueiredo, R. D. Nowak, and S. J. Wright, “Gradient projection for sparse reconstruction: Application to compressed sensing and other inverse problems,” IEEE J. Sel. Top. Signal Process. 1(4), 586–597 (2007). [CrossRef]

17. S. G. Mallat and Z. Zhang, “Matching pursuits with time-frequency dictionaries,” IEEE Trans. Signal Process. 41(12), 3397–3415 (1993). [CrossRef]

18. J. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” IEEE Trans. Inf. Theory 53(12), 4655–4666 (2007). [CrossRef]

19. T. T. Do, L. Gan, N. Nguyen, and T. D. Tran, “Sparsity adaptive matching pursuit algorithm for practical compressed sensing,” In 2008 42nd Asilomar Conference on Signals, Systems and Computers (IEEE, 2008), pp. 581–587. [CrossRef]

20. D. Needell and J. A. Tropp, “CoSaMP: Iterative signal recovery from incomplete and inaccurate samples,” Appl. Comput. Harmon. Anal. 26(3), 301–321 (2009). [CrossRef]

21. J. Tropp, “Greed is good: Algorithmic results for sparse approximation,” IEEE Trans. Inf. Theory 50(10), 2231–2242 (2004). [CrossRef]

22. M. Lustig, D. L. Donoho, J. M. Santos, and J. M. Pauly, “Compressed sensing MRI,” IEEE Signal Process. Mag. 25(2), 72–82 (2008). [CrossRef]

23. D. L. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory 52(4), 1289–1306 (2006). [CrossRef]

24. G. H. Chen, J. Tang, and S. Leng, “Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets,” Med. Phys. 35(2), 660–663 (2008). [CrossRef] [PubMed]

25. R. Tur, Y. C. Eldar, and Z. Friedman, “Innovation rate sampling of pulse streams with application to ultrasound imaging,” IEEE Trans. Signal Process. 59(4), 1827–1842 (2011). [CrossRef]

26. J. Liu, Q. He, and J. Luo, “Compressed sensing for high frame rate, high resolution and high contrast ultrasound imaging,” In 2015 37th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC) (IEEE, 2015), pp. 1552–1555.

27. J. Provost and F. Lesage, “The application of compressed sensing for photo-acoustic tomography,” IEEE Trans. Med. Imaging 28(4), 585–594 (2009). [CrossRef] [PubMed]

28. F. Bunea, “Honest variable selection in linear and logistic regression models via l1 and l 1 + l 2 penalization,” Electron. J. Stat. 2(0), 1153–1194 (2008). [CrossRef]

29. S. M. Kakade, O. Shamir, K. Sridharan, and A. Tewari, “Learning exponential families in high-dimensions: Strong convexity and sparsity,” arXiv preprint arXiv: 0911.0054 (2009).

30. S. Negahban, B. Yu, M. J. Wainwright, and P. K. Ravikumar, “A unified framework for high-dimensional analysis of M -estimators with decomposable regularizers,” In Proceeding of Advances in Neural Information Processing Systems, Y. Bengio, D. Schuurmans, J. D. Lafferty, C. K. I. Williams and A. Culotta, ed. (MIT, 2009), pp. 1348–1356.

31. A. Agarwal, S. Negahban, and M. J. Wainwright, “Fast global convergence rates of gradient methods for high-dimensional statistical recovery,” In Proceeding of Advances in Neural Information Processing Systems, J. D. Lafferty, C. K. I. Williams, J. Shawe-Taylor, R. S. Zemel and A. Culotta, ed. (MIT, 2010), pp. 37–45.

32. T. Blumensath and M. E. Davies, “Gradient pursuit for non-linear sparse signal modelling,” In 2008 16th European Signal Processing Conference (IEEE, 2008), pp. 1–5.

33. F. Dupé, “Greed is Fine: on Finding Sparse Zeros of Hilbert Operators,” in Proceedings of the 31st International Conference on Machine Learning, E. P. Xing and T. Jebara, ed. (Microtome Publ, 2015), Vol. 37.

34. T. Blumensath, “Compressed sensing with nonlinear observations and related nonlinear optimization problems,” IEEE Trans. Inf. Theory 59(6), 3466–3474 (2013). [CrossRef]

35. X. Yuan, P. Li, and T. Zhang, “Gradient hard thresholding pursuit for sparsity-constrained optimization,” arXiv preprint arXiv:1311.5750 (2013).

36. S. Bahmani, B. Raj, and P. T. Boufounos, “Greedy sparsity-constrained optimization,” J. Mach. Learn. Res. 14(1), 807–841 (2013).

37. A. J. Dobson and A. Barnett, An Introduction to Generalized Linear Models (CRC press, 2008).

38. S. Negahban, P. Ravikumar, M. J. Wainwright, and B. Yu, “A unified framework for high-dimensional analysis of M-estimators with decomposable regularizers,” Manuscript, University of California, Berkeley, Dept. of Statistics and EECS (2011).

39. B. Zhang, X. Cao, F. Liu, X. Liu, X. Wang, and J. Bai, “Early-photon fluorescence tomography of a heterogeneous mouse model with the telegraph equation,” Appl. Opt. 50(28), 5397–5407 (2011). [CrossRef] [PubMed]

40. S. R. Arridge, M. Schweiger, M. Hiraoka, and D. T. Delpy, “A finite element approach for modeling photon transport in tissue,” Med. Phys. 20(2), 299–309 (1993). [CrossRef] [PubMed]

41. F. Gao, H. Zhao, L. Zhang, Y. Tanikawa, A. Marjono, and Y. Yamada, “A self-normalized, full time-resolved method for fluorescence diffuse optical tomography,” Opt. Express 16(17), 13104–13121 (2008). [CrossRef] [PubMed]

42. D. Zhu, Y. Zhao, R. Baikejiang, Z. Yuan, and C. Li, “Comparison of regularization methods in fluorescence molecular tomography,” Comparison of regularization methods in fluorescence molecular tomography, in Photonics (Multidisciplinary Digital Publishing Institute, 2014), pp. 95–109.

43. J. C. Baritaux, K. Hassler, M. Bucher, S. Sanyal, and M. Unser, “Sparsity-driven reconstruction for FDOT with anatomical priors,” IEEE Trans. Med. Imaging 30(5), 1143–1153 (2011). [CrossRef] [PubMed]

44. H. Lee, M. Y. Berezin, M. Henary, L. Strekowski, and S. Achilefu, “Fluorescence lifetime properties of near-infrared cyanine dyes in relation to their structures,” J. Photochem. Photobiol. Chem. 200(2-3), 438–444 (2008). [CrossRef] [PubMed]

45. C. A. Mela, C. Patterson, W. K. Thompson, F. Papay, and Y. Liu, “Stereoscopic Integrated Imaging Goggles for Multimodal Intraoperative Image Guidance,” PLoS One 10(11), e0141956 (2015). [CrossRef] [PubMed]

46. G. Alexandrakis, F. R. Rannou, and A. F. Chatziioannou, “Tomographic bioluminescence imaging by use of a combined optical-PET (OPET) system: a computer simulation feasibility study,” Phys. Med. Biol. 50(17), 4225–4241 (2005). [CrossRef] [PubMed]

47. J. Shi, F. Liu, H. Pu, S. Zuo, J. Luo, and J. Bai, “An adaptive support driven reweighted L1-regularization algorithm for fluorescence molecular tomography,” Biomed. Opt. Express 5(11), 4039–4052 (2014). [CrossRef] [PubMed]

Experiment	EED (mm)	Reconstruction method	Evaluation metrics for $I_{τ}$		Evaluation metrics for $τ$
Experiment	EED (mm)	Reconstruction method	MAE ( $10^{- 9} s^{- 1}$ )	CNR	RE_max	Separable
1	9	ANOMP	0.06	4.74	4.90%	Yes
1	9	L1PSD	0.18	2.68	6.77%	No
2	5	ANOMP	0.04	5.52	9.10%	Yes
2	5	L1PSD	0.14	2.72	7.87%	No
3	2	ANOMP	0.02	7.37	4.38%	Yes
3	2	L1PSD	0.08	3.69	27.48%	No

Nonlinear greedy sparsity-constrained algorithm for direct reconstruction of fluorescence molecular lifetime tomography

Abstract

1. Introduction

2. Methods

2.1 Forward model for light propagation

2.2 FMLT reconstruction based on ANOMP algorithm

Parameters setting

ANOMP

2.3 Quantitative metrics for FMLT reconstruction

3. Experiments and results

3.1 Phantom experiments

3.2 In vivo mouse experiments

4. Discussions

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (5)

Tables (1)

Equations (19)

Biomedical Optics Express