Alternating minimization of the negative Poisson likelihood function for the global analysis of fluorescence lifetime imaging microscopy data

Jiyeong Seok; Jeongtae Kim

doi:10.1364/OE.22.024977

1. Introduction

In order to successfully use fluorescence lifetime imaging microscopy (FLIM) to study the molecular environments of a cell, high performance lifetime estimations of decaying fluorescence signals is vital because the lifetime is an indicator of environments [1, 2]. The performance of the lifetime estimation may improve as the photon count increases owing to the signal to noise ratio (SNR) of the decaying fluorescence signals being higher when the photon count is larger [3]. Therefore, one may attempt to increase the exposure time for acquiring high photon count signals, thereby improving the performance of the lifetime estimation. However, the attempt often fails due to photo bleaching and limitations of the acquisition time [3].

Since poor SNR signals often result in noisy lifetime images, several investigations attempted to reduce the effect of the noise in estimated lifetime images using methods such as the penalized maximum likelihood method [2], the wavelet based denoising method [4] and the total variation based denoising method [4, 5]. In addition, denoising methods such as the non-local PCA based method [6] and the sparsity based denoising method using dictionary leaning [7] have possibility to be applied for denoising fluorescence lifetime images. Besides that, the compressive sensing strategy, which was applied to reduce the acquisition time for fluorescence microscopy [8], might be applied for improving the performance of lifetime estimation.

Another very useful method that can improve the performance of the lifetime estimation without increasing exposure time is the global analysis method [3, 9, 10]. The main idea of the global analysis method is that simultaneous using of all the data from the locations within the area, which may significantly improve the performance of the lifetime estimation. This can be accomplished by fitting an exponentially decaying model with the same lifetime but different concentrations at various locations to the measured fluorescence signals from the area. Several studies reported that the global analysis was able to dramatically improve the performance of the lifetime estimation [3, 9, 10].

One of the most serious problems of the global analysis method is the difficulty of optimization since it requires a simultaneous optimization of large amount of variables [10]. For example, if one fits an exponentially decaying model with J lifetimes to the data from I pixels, the global analysis requires the simultaneous optimization of J lifetime and I × J amplitude variables with the non-convex objective function. As the number of variables increases, the computational requirement for the optimization can become more demanding. To improve the speed of the optimization, a previous investigation applied an image segmentation method for a good initial guess [10]. However, even with a good initial guess, optimization may be intractable using the general optimization methods that require the inversion of the Hessian matrix because the computational requirement for the inversion increases almost cubically with the number of variables. In addition, the previous investigation used the χ²-based method, which has been shown to perform worse than the maximum likelihood estimation (MLE) method, which is based on the maximization of the Poisson likelihood function [11]. Another investigation on the rapid global analysis also used the χ²-based objective function to utilize the special structure of the sum of square errors in the objective function [12]. However, the method may not be applied for maximizing the Poisson likelihood function.

In this investigation, to maximize the Poisson likelihood function, we propose an efficient optimization strategy that is based on the alternating minimization method and the optimization transfer approach [13–15]. The alternating optimization strategy has been used effectively for optimization problems that become easier if some variables are fixed [16, 17]. The strategy was also successfully applied for the restoration of Poissonian image [18]. Using this strategy, we alternatively optimize the lifetime and the amplitude variables. We apply the Gauss-Newton method for the optimization of the lifetime variables and the optimization transfer approach for the optimization of the amplitude variables. We note that there exists a previous investigation on the expectation maximization (EM) based estimation of lifetime parameters [19], which has the similar spirit to the optimization transfer approach. However, the EM based method in the previous investigation is limited for estimating parameters in single exponential function. Moreover, it is difficult to apply the method for global analysis [19]. In comparison, the proposed method rapidly fits multi-exponentially decaying functions to acquired signals using global analysis. In simulation studies, we show that the proposed method is able to maximize the Poisson likelihood function much faster than the conventional simultaneous optimization method.

2. Theory

2.1. Problem formulation

We model the number of photon counts detected at a particular location x_i at the time instance t_k using the Poisson random variable, y_i(t_k), as follows [1]:

y_{i} (t_{k}) ~ Poisson {λ_{i} (t_{k}; τ, A_{i})}, k = 0, 1, \dots, K - 1 .

The mean of the Poisson random variable is defined as

λ_{i} (t_{k}; τ, A_{i}) = \sum_{j = 0}^{J - 1} a_{i j} g (t_{k}; τ_{j}),

where a_ij is the amplitude of the j-th decaying function g(t; τ_j), which is defined by the convolution of an exponentially decaying function with the j-th decay constant τ_j and the impulse response function (IRF) of the measurement system, h(t), as follows:

g (t; τ_{j}) = e^{- \frac{t}{τ_{j}}} * h (t) .

We denote the collection of decay constants τ = [τ_j|j = 0,,J − 1] and the collection of amplitudes A_i = [a_ij|j = 0,,J − 1] by the lifetime parameter vector and the concentration parameter vector at x_i, respectively. Note that the lifetime parameter vector τ is the same for every pixel location, while the concentration parameter vector is allowed to have different values at different locations. The log-likelihood function of the entire parameter vector τ and A = {A_i|i = 0,,I − 1} from the entire measurement data is defined as

L (τ, A; y) = \sum_{k = 0}^{K - 1} \sum_{i = 0}^{I - 1} - λ_{i} (t_{k}; τ, A_{i}) + y_{i} (t_{k}) log λ_{i} (t_{k}; τ, A_{i}) - log y_{i} (t_{k})!,

where y_i = [y_i(t_k)|k = 0,,K − 1] is the measured data vector at x_i and y = [y_i|i = 0,,I − 1] is the collection of the measurements from all locations. The maximum likelihood estimation (MLE) of (τ, A) from the measured data y = [y_i|i = 0,,I − 1] is given by

(\hat{τ}, \hat{A}) = \underset{τ, A > 0}{argmin} - L (τ, A; y) .

Note that the minimization problem defined in Eq. (5) is a coupled optimization problem because τ is a common parameter vector for every pixel location. Therefore, solving the optimization problem requires the simultaneous optimization of a large sized variable. One may attempt to solve the optimization problem using a software package such as the Optimization toolbox in MATLAB^® (Mathworks, USA), as in a previous investigation [10]. However, the computation time can be demanding as the size of variables increases because the inversion of the large size Hessian matrix is often involved. Note that the size of the variable increases as more data points are simultaneously analyzed.

We note that there exist investigations based on different noise model from the one adopted in this research. For example, several investigations used the space variant Gaussian noise model [5, 20] based on the approximation of the Poisson distribution to the Gaussian distribution. However, it was shown that the performance of lifetime estimation using such approximation was worse than that of the Poisson noise based method [2]. There also exist investigations based on the sum of the Poisson and the Gaussian noise model [21, 22]. This model might be more accurate than ours since it considers both random photon arrivals and electric circuit noise. However, because the Poisson noise is dominant over the Gaussian noise in low photon signals (such as signals that require global analysis), we mainly focus on the Poisson noise in this investigation. We note that our optimization strategy in this investigation can be easily extended for the sum of the Poisson and the Gaussian noise model as it was done previously in medical image reconstruction [14].

2.2. Proposed method

To design an efficient method to determine the minimizer of the negative Poisson log-likelihood function defined in Eq. (4), reducing the variable size during the optimization is desirable to avoid inverting the large size Hessian matrix. In addition, if the lifetime parameters are fixed, then the minimization of the negative Poisson likelihood function with respect to the concentration parameters is a convex optimization problem, which may be handled more easily than the non-convex optimization problem. These two facts suggest the following alternating optimization method that iteratively solves two optimization problems, which are the minimizations of the negative log-likelihood function with respect to the lifetime parameters (concentration parameters) while fixing the concentration parameters (lifetime parameters) using the most recently determined values, as follows:

τ^{n + 1} = \underset{τ > 0}{argmin} - L (τ; A^{n}, y),

A^{n + 1} = \underset{A > 0}{argmin} - L (A; τ^{n + 1}, y) .

Because the dimensionality of the lifetime parameter vector is usually small, e.g., J for an exponential decay model with J different lifetimes, the computational complexity to determine the inverse of the Hessian matrix for the optimization problem defined in Eq. (6) is not high. We solved the optimization problem using the Gauss-Newton method defined as follows:

τ^{n + 1} = τ^{n} + H^{- 1} (τ^{n}; A, y) \nabla L (τ^{n}; A^{n}, y),

where we approximated the J × J Hessian matrix H(τⁿ; A, y) = [h_pq] as

h_{p q} = \sum_{k = 0}^{K - 1} \sum_{i = 0}^{I - 1} (\frac{y_{i} (t_{k})}{{(\sum_{j} a_{i j} g (t_{k}; τ_{j}^{n}))}^{2}}) (\frac{a_{i j} t_{k}}{{τ_{p}}^{2}} g (t_{k}; {τ_{p}}^{n})) (\frac{a_{i j} t_{k}}{{τ_{q}}^{2}} g (t_{k}; {τ_{q}}^{n})) .

We ignored the second order derivative terms when calculating the Hessian matrix because including these terms may destabilize the optimization [11, 23]. Unlike the optimization problem defined in Eq. (6), the optimization of the concentrations defined in Eq. (7) is a convex optimization problem with large amount of variables. Because the size of the concentration parameters is usually much larger than the size of lifetime parameters, it may not be effective to apply a general optimization method that requires the inversion of the Hessian matrix. We solve the optimization problem by applying an optimization transfer strategy, which is based on the fact that the sequence of minimizers of surrogate functions Q(θ; θ^m), in which θ^m denotes the estimated value of θ at the m-th iteration, converges to the minimizer of the original objective function Φ(θ) if the surrogate functions satisfy the following two conditions [14, 15]:

Q (θ; θ^{m}) \geq Φ (θ),

Q (θ^{m}; θ^{m}) = Φ (θ^{m}) .

If the minimizer of the surrogate function can be easily determined, then minimizing the sequences of surrogate functions can be more effective than the direct minimization of the original objective function. The optimization transfer approach has been successfully applied for image restoration and reconstruction problems, which require simultaneous optimizations of large sized variables for a convex objective function [14,24]. We design an additively separable surrogate function for the objective function in Eq. (7) to determine

a_{i j}^{n + 1}

. We construct a separable surrogate function so that the minimizer of the surrogate function can be independently determined. Furthermore, we design the surrogate function such that the minimizer of the surrogate function can be determined in a closed form without using an iterative optimization. To construct such a separable surrogate function, we apply De Pierros’s method using the Jansen’s inequality [24]. Note that the log-likelihood function in Eq. (7) can be represented as

- L (A; τ^{n + 1}, y) = \sum_{k = 0}^{K - 1} \sum_{i = 0}^{I - 1} {\sum_{j = 0}^{J - 1} a_{i j} g (t_{k}; τ_{j}^{n + 1}) - y_{i} (t_{k}) log (\sum_{j = 0}^{J - 1} a_{i j} g (t_{k}; τ_{j}^{n + 1}))} .

First, we note the following relationship,

\sum_{j} a_{i j} g (t_{k}; τ_{j}^{n + 1}) = \sum_{j} (\frac{a_{i j}^{(n, m)} g (t_{k}; τ_{j}^{n + 1})}{\sum_{j} a_{i j}^{(n, m)} g (t_{k}; τ_{j}^{n + 1})}) (\frac{a_{i j}}{a_{i j}^{(n, m)}} \sum_{j} a_{i j}^{(n, m)} g (t_{k}; τ_{j}^{n + 1})),

where

a_{i j}^{(n, m)}

is the estimated value of a_ij during the m-th sub-iteration to determine

a_{i j}^{(n + 1)}

. Substituting the argument of the logarithm function in Eq. (12) into Eq. (13) and applying the Jansen’s inequality yields the following inequality:

\begin{array}{l} - log (\sum_{j} (\frac{a_{i j}^{(n, m)} g (t_{k}; τ_{j}^{n + 1})}{\sum_{j} a_{i j}^{(n, m)} g (t_{k}; τ_{j}^{n + 1})}) (\frac{a_{i j}}{a_{i j}^{(n, m)}} \sum_{j} a_{i j}^{(n, m)} g (t_{k}; τ_{j}^{n + 1}))) \\ \leq - \sum_{j} (\frac{a_{i j}^{(n, m)} g (t_{k}; τ_{j}^{n + 1})}{\sum_{j} a_{i j}^{(n, m)} g (t_{k}; τ_{j}^{n + 1})}) log (\frac{a_{i j}}{a_{i j}^{(n, m)}} \sum_{j} a_{i j}^{(n, m)} g (t_{k}; τ_{j}^{n + 1})) . \end{array}

Furthermore, the equality in Eq. (14) holds if

a_{i j} = a_{i j}^{(n, m)}

. Therefore, the following function can be a valid surrogate function:

Q (A; A^{(n, m)}) = \sum_{i, j, k} a_{i j} g (t_{k}; τ_{j}^{n + 1}) - (\frac{y_{i} (t_{k}) a_{i j}^{(n, m)} g (t_{k}; τ_{j}^{n + 1})}{\sum_{j} a_{i j}^{(n, m)} g (t_{k}; τ_{j}^{n + 1})}) log (\frac{a_{i j}}{a_{i j}^{(n, m)}} \sum_{j} a_{i j}^{(n, m)} g (t_{k}; τ_{j}^{n + 1})) .

The next estimate can be determined using the following derivative with respect to a_ij:

\frac{\partial Q (A; A^{(n, m)})}{\partial a_{i j}} = \sum_{k} g (t_{k}; τ_{j}^{n + 1}) - \sum_{k} y_{i} (t_{k}) (\frac{a_{i j}^{(n, m)} g (t_{k}; τ_{j}^{n + 1})}{\sum_{j} a_{i j}^{(n, m)} g (t_{k}; τ_{j}^{n + 1})}) \frac{1}{a_{i j}} .

By equating the gradient with zero, the next estimate is determined as follows:

a_{i j}^{(n, m + 1)} = \frac{a_{i j}^{(n, m)}}{\sum_{k} g (t_{k}; τ_{j}^{n + 1})} \sum_{k} (\frac{y_{i} (t_{k}) g (t_{k}; τ_{j}^{n + 1})}{\sum_{j} a_{i j}^{n} g (t_{k}; τ_{j}^{n + 1})}) .

One may perform this sub-iteration M times to define

a_{i j}^{n + 1}

. In other words,

a_{i j}^{n + 1} = a_{i j}^{(n, M)} .

Note that as M approaches infinity,

a_{i j}^{n, M}

converges to the minimizer of the objective function defined in Eq. (7). However, in practice, a few sub-iterations that increase the log-likelihood function may be sufficient.

2.3. Initial estimates

It is important to select initial parameters judiciously to avoid possible local minima and achieve fast convergence. We determine the initial lifetime parameters by the results of a time invariant fit based on the sum of observed signals, which was shown to be more accurate than pixel-wise independent lifetime estimation methods [10]. We initialize the concentration parameters at each pixel location using the maximum value of observed signal divided by the number of exponential functions in the signal model.

3. Results

3.1. Simulation 1

To compare the performance of the proposed method with the simultaneous optimization strategy, we conducted simulation studies using Poisson random variables y_i(t_k), the means of which are double-exponentially decaying values defined as

λ_{i} (t_{k}; τ_{1}, τ_{2}, c_{i 1}, A) = A * {c_{i 1} g (t_{k}; τ_{1}) + (1 - c_{i 1}) g (t_{k}; τ_{2})},

where A is an intensity value and c_i1 is the intensity ratio of the lifetime τ₁ at the i-th pixel location. Note that a similar model to the one defined in Eq. (19) was used in a previous investigation [10]. In the model defined in Eq. (19), the product of A and c_ij is equivalent to a_ij in the model defined in Eq. (3). We generated a 32 × 32 size synthetic image based on Eq. (19) with c_i1 varying from 0 to 1. Figure 1(a) shows the true c_i1 for every pixel location. We set the true values to τ₁ = 2 ns and τ₂ = 4 ns and used the IRF of the Gaussian function with σ = 1 and 5 nonzero points. The average total photon counts at each pixel location was 2938.

Fig. 1 Intensity parameter images from different optimization methods: (a) true intensity parameter image; (b) estimated intensity parameter image using pixel-wise independent optimization without global analysis; (c) estimated intensity parameter image using the trust region reflective simultaneous optimization method; (d) estimated intensity parameter image using the proposed method

Download Full Size | PDF

First, we attempted to recover the lifetime and amplitude parameters without a global analysis. This was done by maximizing the Poisson log-likelihood function at each pixel independently. Next, we applied the global analysis method to improve the performance of the lifetime and amplitude parameters estimation. To determine the minimizer of the negative Poisson likelihood function defined in Eq. (4), we simultaneously optimized all parameters using the f mincon function in the Optimization toolbox in MATLAB^®, which is based on the trust region reflective method [25–27]. We further applied the proposed method to solve the optimization problem defined in Eq. (4) with an improved speed. We ran these optimizations on a workstation equipped with 12 core Intel Xeon E5 processor (2.7 GHz) and 64 GB of memory.

Figure 1(b) through Fig. 1(d) show the estimated c_i1 images using the pixel-wise independent Poisson likelihood function maximization method, global analysis using simultaneous optimization method, and global analysis using the proposed optimization method, respectively. As shown in Fig. 1(b), the independent Poisson likelihood maximization at each pixel yielded a noisy image due to the poor SNR. The global analysis method was able to improve the performance of the estimation as shown in Fig. 1(c) and Fig. 1(d), which are closer to the true image shown in Fig. 1(a). Both the simultaneous optimization method and proposed method generated nearly identical results, implying that the two methods found the same minimum point during the optimization. However, the proposed method determined the minimum point significantly faster than the conventional simultaneous optimization method. Figure 2(a) reports the convergence of the negative Poisson likelihood function for the global analysis using the simultaneous optimization and proposed alternating optimization method as the function of central processing unit (CPU) time. As shown in Fig. 2(a), the proposed alternating optimization method decreases the Poisson likelihood function much faster than the simultaneous optimization method. The CPU time required to reach the negative Poisson likelihood function value of −5.837 × 10⁶, the almost converged value, using the proposed method was approximately 1 seconds, while that required by the conventional method was approximately 1050 seconds. Note that the proposed method was able to decrease the objective function monotonically without using a line search algorithm. Figure 2(b) through Fig. 2(e) show the change of parameters (τ₁, τ₂, a_i1, and a_i2) during the optimization. Since the amplitude value at each pixel location is different, we report parameter values at one pixel location: the (8, 16) location. The true c_i1 value of the pixel location was 0.2424. Therefore, the amplitude of the first exponential function was a_i1 = 50 × c_i1 = 12.12 and that for the second exponential function was a_i2 = 50 × (1 − c_i1) = 37.88. As shown in Fig. 2, the proposed method converged significantly faster than the conventional simultaneous optimization. Note that even in the time interval where the change of the negative Poisson likelihood function is very small, after approximately 500 s CPU time, the parameters were still changing in the conventional method. As shown in Fig. 2, both the proposed method and simultaneous optimization method converged to nearly the same parameter values (τ₁ = 2.0026 ns, τ₂ = 4.0072 ns, A₁ = 17.49, and A₂ = 35.42). However, as mentioned before, the convergence of the proposed alternating optimization method was much faster (by at least 100 times) than the simultaneous optimization method.

Fig. 2 Likelihood value and estimated parameters as functions of CPU time for different optimization methods: (a) likelihood; (b) τ₁; (c) τ₂; (d) A₁ of (8, 16) location; (e) A₂ of (8, 16) location;

Download Full Size | PDF

3.2. Simulation 2

We conducted another simulation study using Poisson random variables y_i(t_k), whose means are tri-exponentially decaying signals defined as follows:

λ_{i} (t_{k}; τ_{1}, τ_{2}, A_{1}, A_{2}, A_{3}) = A_{1 i} g (t_{k}; τ_{1}) + A_{2 i} g (t_{k}; τ_{2}) + A_{3 i} g (t_{k}; τ_{2}),

where the decay constants are the same (τ₁=1ns, τ₂=2ns, τ₃=4ns) for every pixel but associated amplitudes (A_1i, A_2i, A_3i) may be different at different pixel location. The number of data points at each pixel location was 256. Figs. 3(a)–3(c) shows the images of true amplitude 1, 2 and 3, respectively. We attempted to estimate the decay constants and the amplitudes by fitting bi-exponentially decaying functions to noisy tri-exponentially decaying signals. Note that this sort of modeling error may happen in real experiments since the number of exponential functions which are present in acquired signals may not always be available.

Fig. 3 True amplitude images in simulation 2: (a) amplitude 1; (b) amplitude 2; (c) amplitude 3;

Download Full Size | PDF

Figs. 4(a) and 4(b) show the convergence of the two decay constants during optimization. As one can see in the figures, the convergence of the proposed method is much faster than that of the conventional simultaneous optimization method while the two methods determine almost identical decay constants, τ₁=1.3ns and τ₂=3.76ns. Note that the presence of un-modeled exponential function results in inaccurate estimation of the decay constants. Figs. 4(c) and 4(d) show estimated amplitude images using the global analysis while Figs. 4(e) and 4(f) show those without using the global analysis strategy. As in the previous simulation, the amplitude images determined using the proposed method and the simultaneous global analysis were almost identical. However, the pixel-wise independent estimation yielded very noisy image. In this case, estimated decay constants were 2.76 ns and 29.36ms, which were very inaccurate. We note that the purpose of this simulation is to show that the proposed method converges much faster than the global method with simultaneous optimization even if there exist errors in signal modeling.

Fig. 4 Estimated decay parameters as functions of CPU time for different optimization methods: (a) τ₁; (b) τ₂; (c) amplitude 1 using the global analysis; (d) amplitude 2 using the global analysis; (e) amplitude 1 using the pixel-wise independent estimation; (f) amplitude 2 using the pixel-wise independent estimation;

Download Full Size | PDF

4. Conclusion

We investigated a fast optimization method for the global analysis for FLIM. Unlike the existing simultaneous optimization of concentration and lifetime parameters, the proposed method is based on the alternating optimization strategy, which makes the optimization of concentration parameters a convex optimization problem. We effectively determined the minimizer of the convex objective function using the optimization transfer strategy based on the convex inequality, and we solved a non-convex optimization problem to estimate the lifetime parameters using the Gauss-Newton method. In the simulation studies, the proposed method determined nearly the same parameters as the ones found by the conventional simultaneous optimization method, within a significantly reduced time. As a result, we believe that the proposed optimization method is very useful for the global analysis of FLIM.

Acknowledgments

This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology ( NRF-2012R1A1A2009370) and Ewha Global Top 5 Grant 2011 of Ewha Womans University.

References and links

1. J. R. Lakowicz, Principles of Fluorescence Spectroscopy (Kluwer Academic/Plenum, 1999). [CrossRef]

2. J. Kim, J. Seok, H. Lee, and M. Lee, “Penalized maximum likelihood estimation of lifetime and amplitude images from multi-exponentially decaying fluorescence signals,” Opt. Express 21, 20240–20253 (2013). [CrossRef] [PubMed]

3. P. J. Verveer, A. Squire, and P. I. Bastiaens, “Global analysis of fluorescence lifetime imaging microscopy data,” Biophys. J. 78, 2127–2137 (2000). [CrossRef] [PubMed]

4. C. W. Chang and M.-A. Mycek, “Total variation versus wavelet-based methods for image denoising in fluorescence lifetime imaging microscopy,” J. Biophotonics 5, 449–457 (2012). [CrossRef] [PubMed]

5. C. W. Chang and M.-A. Mycek, “Enhancing precision in time-domain fluorescence lifetime imaging,” J. Biomed. Opt. 15, 056013 (2010). [CrossRef] [PubMed]

6. J. Salmon, C.-A. Deledalle, R. Willett, and Z. Harmany, “Poisson noise reduction with non-local PCA,” in proceedings of IEEE International Conference on Acoustics, Speech and Signal Processing (IEEE, 2012), pp. 1109–1112.

7. R. Giryes and M. Elad, “Sparsity based poisson denoising,” in proceedings of IEEE Convention of Electrical Electronics Engineers in Israel (IEEE, 2012), pp. 1–5.

8. V. Studer, J. Bobin, M. Chahid, H. S. Mousavi, E. Candes, and M. Dahan, “Compressive fluorescence microscopy for biological and hyperspectral imaging,” Proceedings of the National Academy of Sciences 109, E1679–E1687 (2012). [CrossRef]

9. H. E. Grecco, P. Roda-Navarro, and P. J. Verveer, “Global analysis of time correlated single photon counting FRER-FLIM data,” Opt. Express 17, 6493–6508 (2009). [CrossRef] [PubMed]

10. S. Pelet, M. J. R. Previte, L. H. Laiho, and P. T. C. So, “A fast global fitting algorithm for fluorescence lifetime imaging microscopy based on image segmentation,” Biophys. J. 87, 2807–2817 (2004). [CrossRef] [PubMed]

11. J. Kim and J. Seok, “Statistical properties of amplitude and decay parameter estimators for fluorescence lifetime imaging,” Opt. Express 21, 6061–6075 (2013). [CrossRef] [PubMed]

12. S. C. Warren, A. Margineanu, D. Alibhai, D. J. Kelly, C. Talbot, Y. Alexandrov, I. Munro, M. Katan, C. Dunsby, and P. M. W. French, “Rapid global fitting of large fluorescence lifetime imaging microscopy datasets,” PLoS ONE 8, e70687 (2013). [CrossRef] [PubMed]

13. J. Fessler, E. Ficaro, N. Clinthorne, and K. Lange, “Grouped-coordinate ascent algorithms for penalized-likelihood transmission image reconstruction,” IEEE Trans. Med. Imag. 16, 166–175 (1997). [CrossRef]

14. J. Fessler, “Image reconstruction: Algorithms and analysis,” Online preprint of book in preparation.

15. K. Lange, Optimization, Springer Texts in Statistics (Springer, 2004). [CrossRef]

16. U. Niesen, D. Shah, and G. W. Wornell, “Adaptive alternating minimization algorithms,” IEEE Trans. Inf. Theory 55, 1423–1429 (2009). [CrossRef]

17. Y. Wang, J. Yang, W. Yin, and Y. Zhang, “A new alternating minimization algorithm for total variation image reconstruction,” SIAM J. Imaging Sci.248–272 (2008). [CrossRef]

18. M. Figueiredo and J. Bioucas-Dias, “Restoration of poissonian images using alternating direction optimization,” IEEE Trans. Image Process. 19, 3133–3145 (2010). [CrossRef] [PubMed]

19. A. Jezierska, C. Chaux, J.-C. Pesquet, H. Talbot, and G. Engler, “An EM approach for time-variant poisson-gaussian model parameter estimation,” IEEE Trans. Signal Process. 62, 17–30 (2014). [CrossRef]

20. C. W. Chang and M.-A. Mycek, “Precise fluorophore lifetime mapping in live-cell, multi-photon excitation microscopy,” Opt. Express 18, 8688–8696 (2010). [CrossRef] [PubMed]

21. P. Roudot, C. Kervrann, F. Waharte, and J. Boulanger, “Lifetime map reconstruction in frequency-domain fluorescence lifetime imaging microscopy,” in proceedings of IEEE International Conference on Image Processing, (IEEE, 2012), pp. 2537–2540.

22. P. Roudot, C. Kervrann, J. Boulanger, and F. Waharte, “Noise modeling for intensified camera in fluorescence imaging: Application to image denoising,” in proceedings of IEEE International Symposium on Biomedical Imaging (IEEE, 2013), pp. 600–603.

23. W. H. Press, Numerical recipes : the art of scientific computing, 3rd ed. (Cambridge University, 2007).

24. A. De Pierro, “A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography,” IEEE Trans. Med. Imag. 14, 132–137 (1995). [CrossRef]

25. J. Moré and D. Sorensen, “Computing a trust region step,” SIAM J. Sci. Comput. 4, 553–572 (1983). [CrossRef]

26. M. A. Branch, T. F. Coleman, and Y. Li, “A subspace, interior, and conjugate gradient method for large-scale bound-constrained minimization problems,” SIAM J. Sci. Comput. 21, 1–23 (1999). [CrossRef]

27. R. Byrd, R. Schnabel, and G. Shultz, “Approximate solution of the trust region problem by minimization over two-dimensional subspaces,” Math. Program. 40, 247–263 (1988). [CrossRef]

Alternating minimization of the negative Poisson likelihood function for the global analysis of fluorescence lifetime imaging microscopy data

Abstract

1. Introduction

2. Theory

2.1. Problem formulation

2.2. Proposed method

2.3. Initial estimates

3. Results

3.1. Simulation 1

3.2. Simulation 2

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (4)

Equations (20)

Optics Express