Blind and reference-free fluorescence lifetime estimation via consumer time-of-flight sensors

Fluorescence lifetime imaging (FLI) is a popular method for extracting useful information that is otherwise unavail-able from a conventional intensity image. Usually, however, it requires expensive equipment, is often limited to either distinctly frequency- or time-domain modalities, and demands calibration measurements and precise knowledge of the illumination signal. Here, we present a generalized time-based, cost-effective method for estimating lifetimes by repur-posing a consumer-grade time-of-flight sensor. By developing mathematical theory that unifies time- and frequency-domain approaches, we can interpret a time-based signal as a combination of multiple frequency measurements. We show that we can estimate lifetimes without knowledge of the illumination signal and without any calibration. We experimentally demonstrate this blind, reference-free method using a quantum dot solution and discuss the method ’ s implementation in FLI applications.


INTRODUCTION
Fluorescence lifetime imaging (FLI) is a significant research area that spans many engineering applications. Knowledge of a sample's fluorescence lifetime allows, for example, DNA sequencing [1], tumor detection [2,3], fluorescence tomography [4,5], in vivo imaging [6], and high-resolution microscopy [7]. Typically, FLI is categorized into two complementary modes [4,8]. In timedomain FLI [7,9] (TD-FLI), an impulse-like excitation pulse probes the fluorescent sample, and the time-resolved reflection is used to calculate lifetimes. In frequency-domain FLI [10][11][12] (FD-FLI), the sample is excited with (sinusoidal) intensity modulated light, and the measured phase shift of the reflected signal at the same modulation frequency encodes the lifetime. FD-FLI is theoretically appealing in that phase measurements at one given modulation frequency suffice to resolve the sample lifetime. From a practical standpoint, model mismatch [10] and sample contamination due to multiple lifetimes [13] often limit the accuracy of single-frequency-based systems. Frequency diversity [13,14] may be used for imaging multiple lifetimes. Since this approach requires sweeping frequencies over a given bandwidth, a frequency diversity method is not suitable for wide-field imaging of dynamic samples. In either case (TD-FLI or FD-FLI), because of the costly equipment, system constraints are often strict, so that precise knowledge of the illumination signal and calibration measurements (to compensate for path length delays) are required [15].
Alternatively, time-of-flight (ToF) sensors, which are far more cost-effective and operate in real time, offer a great deal of flexibility in design. ToF sensors, such as the Microsoft Kinect, are essentially depth/range imaging devices and are at the heart of the entertainment industry. Recently, ToF sensing has found prodigious use in computer graphics, computer vision [16], and computational imaging [16][17][18], with applications in multipath imaging [19][20][21][22], ultrafast imaging [22,23], and imaging through scattering media [17,24].
Here, we demonstrate an extension of ToF sensing to FLI that is both blind and calibration-free. To do so, we generalize (from the distinct time or frequency domains) the active illumination signal, which is common to both modalities, to show that neither calibration measurements nor knowledge of the signal are necessary for estimation of fluorescent lifetimes. This blind, referencefree method is implemented in a consumer-grade ToF sensor to estimate simultaneously the range and lifetime of a CdSe-CdS quantum dot sample. This method suggests a cost-effective alternative to the usual FLI methods. This paper is organized as follows: starting from first principles, we propose a ToF-sensor-based image formation model and discuss its role in depth/range estimation in Section 2. We then discuss the two complementary modes of operation of ToF sensors: time and frequency domains. Section 2.A deals with TD-ToF-sensor-based depth imaging, while Section 2.B discusses frequency-domain depth imaging. In Section 3, we provide a mathematical model for FLI with ToF sensors. We show that consumer ToF sensors can be used to estimate the lifetime of a fluorescent sample. We develop theoretical models for both time-and frequency-domain modes, and our theory is corroborated with experiments for both. These complementary approaches are discussed in Sections 3.A and Section 3.B. Finally, we conclude this work with possible future directions in Section 5.

TOF IMAGE FORMATION MODEL
ToF sensors operate using the lock-in principle [25]: active illumination probes a sample, and the reflected light is cross-correlated electronically to calculate range and amplitude information. In this way, each ToF sensor exposure results in two images: the usual intensity image and a range image in which each pixel relates to the depth of the scene. Usually, the signal is sinusoidal, essentially classifying a ToF sensor as a homodyne detection system, but the illumination signal can be far more general. Both ToF sensors and current FLI technologies acquire scene information in either the time or frequency domains, so the present analysis for ToF sensing is compatible with both.
Regardless of the operating domain or physical implementation, the ToF imaging process can be understood as follows. A scene is illuminated with a time-dependent Δ-periodic probing intensity signal pt, such that pt Δ pt, Δ > 0. Similar to FLI, TD-ToF [21,22] systems use a time-localized pulse pt (not necessarily an impulse); the FD-ToF counterpart uses a modulated probing function pt 1 cosω 0 t, where ω 0 (usually in the megahertz range) is the modulation frequency. Below, however, we do not use a specific form for pt.
This signal interacts with the scene, whose response is characterized by a time-dependent scene response function (SRF) ht; t 0 , where t 0 is a time-domain variable which models the time-variant SRF. This interaction results in the reflected signal rt: rt Z pt 0 ht; t 0 dt 0 : A ToF sensor detects rt and electronically cross correlates it with pt: where pt p−t. The K stored measurements are digital samples of this cross correlation, m k mkT s , k 0; …; K − 1 (T s > 0 is the sampling rate). In many cases of interest, the SRF is shift-invariant, that is (see [26], Sec. 2.2.4, p. 50), so that Eq. (1) represents the convolution operation defined by In the case of shift-invariant SRFs, the measurements simplify to implying linear filtering of cross-correlation function ϕt with the scene response filter h SI . Conventional ToF sensors [16] are designed for range imaging and depth estimation. For simple reflections from an object (with reflection coefficient ρ) that is at a depth d from the sensor, the SRF becomes ht; t 0 ρδt − t 0 − 2d ∕c; (5) where c is the speed of light, and δ is the Dirac distribution. The reflected signal in this case is simply a delayed version of the probing function and the delay is proportional to the depth parameter d. More precisely, in view of the ToF sensor operation, the reflected signal in Eq. (1) reads rt ρpt − t 0 , where the delay t 0 2d ∕c. In this case, the measurements amount to Estimation of the scene parameters fρ; dg by the ToF sensor results in range and intensity images. The estimation process depends on the choice of probing function p, which may be a time-localized pulse or an amplitude-modulated continuous-wave (AMCW) function leading to time-domain [ Fig. 2(a)] and frequency-domain [ Fig. 2(b)] modes of operation, respectively.

A. TD-ToF Imaging
For TD-ToF imaging, the probing function is ideally well localized in time, that is, a Dirac delta distribution. In this case, pt ∼ δt, and the measurements in Eq. (4) simplify to mt h SI t. In practice, maximum length sequences (MLSs) [27] provide an optimal time-localized probing function. Such sequences are a function of the signal length or period Δ.
In view of the depth estimation problem, where ht; t 0 ρδt − t 0 − 2d ∕c, the depth is estimated by the operatioñ More sophisticated methods have been developed for the case of multiple reflections or multipath interference [21,22].
The exact characterization of the true TD-ToF probing function pt is challenging due to the electronics of the sensor and the physical process involved. By using a Fourier series expansion [21] to represent the probing function, The Fourier series coefficients of the cross-correlated signal and the probing signal are related bŷ ϕ n p np n jp n j 2 : Practically, a band-limited approximation, C p;p t, suffices to represent ϕt: n e jnω 0 t ; where we have retained the first N 0 Fourier coefficients of ϕt. Figure 1 shows an experimentally obtained ϕt and its Fourier approximation C p;p , with their Fourier spectra with N 0 30, for an impulsive probe function.

B. FD-ToF Imaging
In FD-ToF imaging, the scene is probed with a continuous wave, sinusoidal probing function with a modulation frequency ω: where p 0 is the modulation amplitude. For pure depth estimation, with the SRF specified in Eq. (5), the reflected signal becomes The ToF lock-in sensor acts as a homodyne detector array: each pixel cross correlates the reflected signal with the probing signal to produce measurements: where, in Eq. (10), we generalize the cross correlation in Eq. (2) for sinusoids. The two quantities of interest, ρ and d , are estimated with the four digital measurements of the measured signal in Eq. (11), that is, m k mπk∕2ωT s , k 0; …; 3. Based on these discrete measurements, we define a complex number z ∈ C: The scene parameters are thus estimated bỹ Note that this technique for extracting phase is identical to conventional phase-shifting holography [28].
To summarize, object range is estimated via time delays in TD-ToF  Table 1 and Table 2, respectively.

Research Article
In most consumer ToF sensors, mt [Eq. (2)] is a set of measurements that is used to estimate the scene parameters fρ; dg [22,25]. The goal of this paper is to show that by using exactly the same set of measurements but with a different SRF, one can recover fluorescent lifetimes in context of the FLI.

FLI WITH TOF SENSORS
The experimental setup for FLI using a ToF sensor is depicted in Fig. 2(e). A fluorescent sample is located in an x-y plane a distance z d from the ToF sensor. Pixels corresponding to a nonfluorescent background object x b ; y b produce only a time delay proportional to 2d ∕c, precisely the case of conventional ToF imaging with the SRF specified in Eq. (5).
The probing function that interacts with the pixels corresponding to a fluorescent sample at location x f ; y f undergoes two transformations. The first transformation is attributed to the same depth contribution, d . The second transformation results from fluorescence: a fraction of the incident light excites the sample, which fluoresces with a characteristic decay time τ. The total SRF at x f ; y f is then where the contributions due to respective components are μ is the emission strength of the fluorescent sample, and Πt is the unit step function. (Note that h Depth and h Sample represent signals at the excitation and emission optical wavelengths, respectively. The former can be eliminated experimentally via, e.g., a dielectric interference filter.) Because Eq. (13) is linear and shiftinvariant, the resulting measurements are We emphasize that our ToF system is completely compatible with FLI. In TD-FLI, pt δt, whereas in FD-FLI, pt 1 cosω 0 t. Conventionally, separate calibration steps provide explicit knowledge of d , which is used for subsequent measurements. However, we make no such assumption. Instead, we simultaneously compute d and τ from our measurements [Eq. (16)]. We need not perform a separate measurement to obtain d explicitly.
A. TD-ToF FLI: Theory and Experiments

Theoretical Modeling
As in TD-ToF, we utilize the same truncated cross-correlation function [Eq. (9)] for FLI. Importantly, explicit knowledge of pt is not required. To show this, note that the eigenfunctions of Eq. (16) are precisely the complex exponentials of Eq. (9). By using the convolution theorem, we have whereĥ n ĥω 0 n andĥω is the Fourier transform of h SI t, From here onward and for all practical purposes, we assume that Eq. (17) is an equality instead of an approximation. Expressingĥω in polar form, we haveĥω jĥωje j∠ĥω , where Note that the phase of the spectrum encodes both the depth and the lifetime parameters. We may write where θ τ ω tan −1 ωτ 2 τ ρ1 ωτ 2 and θ d ω 2d ω∕c: (21) (Because we consider a single lifetime, we set μ 1.) In vector-matrix notation, the discretized system of equations in Eq. (17) can be written as m VDφĥ; where • m is a K × 1 vector of discretized ToF sensor measurements m k mkT s , k ∈ 0; K − 1; • V is a Vandermonde matrix of size K × 2N 0 1 with matrix element V k;n e j2π∕ΔT s nk , n ∈ −N 0 ; N 0 ; • Dφ is a 2N 0 1 × 2N 0 1 diagonal matrix with diagonal entries D n;n φ n ; and •ĥ is 2N 0 1 × 1 vector of the discretized spectrum [Eq. (19)] with entriesĥ2πn∕Δ, n −N 0 ; …; N 0 .
The estimation problem is thus: given K measurements m, estimate parameters d and τ.
Because ϕt is a real, time-symmetric function by construction, the matrix Dφ does not contribute to the phase of vectorĥ in Eq. (22). Hence, we rely on only the measurements m. Thus, p (or ϕ) need not be known, so that we eliminate the calibration requirements that are central to both FLI and ToF imaging [22].
To see this, note that TD-FLI uses [8] p δt ϕ, so the corresponding measurements are Research Article Now because logm TD-FLI t is linear, a direct fit [8] can be used to estimate τ, and, hence, calibration is implicitly avoided by the choice of p δ.
On the other hand, for our proposed TD-ToF method for FLI, the illumination/probing function is composed of N 0 multiplexed frequencies, which yield measurements m TD-ToF-FLI t 17 1 Δ and because C p;p is a band-limited approximation of ϕ [Eq. (9)], a linear fit no longer suffices. However, because of the properties of (the generalized) pt, we solve this more complex inverse problem to estimate fd ; τg without knowing ϕ.

Experimental Verification of TD-ToF FLI
The setup for the TD-ToF FLI method is shown in Fig. 2 with Δ 309.9 ns and T s 7.8120 ps. The sample consists of a CdSe-CdS quantum dot sample prepared by dissolving it in hexane and polymethyl methacrylate (PMMA) onto a glass slide. This sample has a lifetime of τ 32 ns, and it is located at 1.05 m from the sensor. The reflected light is passed through a dielectric interference filter with a cut-off wavelength of 450 nm, leading to ρ 0 [Eq. (21)]. A 160 × 120 pixel PMD 19k-S3 lock-in sensor with custom field programmable gate array (FPGA) programming cross correlates the reflected signal to produce measurements m. The sensor has an 80 MHz modulation bandwidth operating within 90 frames per second [29]. The total cost of our system is $1200.
With N 0 15, we compute ∠ĥ obs ∠V m −θ τ 2πn∕Δ θ d 2πn∕Δ; (26) where V denotes the pseudo-inverse of the matrix V and ∠ĥ obs is the experimentally observed phase. Since we know that the observed phase is related to the theoretical phase in Eq. (20), we estimate the parameters of interest by solving the nonlinear least squares problem arg min where α and β are offset parameters to ensure a solution that is centered at origin. We use the trust-region-based [30] algorithm with the least absolute residual criterion. The nonlinear parameter estimation problem is solved using bounded constraints that are accommodated within the trust-region optimization framework. We require that the minimum and maximum lifetimes and distances are in ranges of 0 to 100 ns and 0 to 10 m, respectively. Per-pixel fluorescent sample measurements are shown in Fig. 3, which shows both the time-domain measurement [ Fig. 3(a)] and the raw phase measurement [ Fig. 3(b)]. We also show the fitted phase measurement in Fig. 3(b). We assign a confidence level to each pixel to estimate the signal-to-noise ratio (SNR). We use ToF measurements from 14 different pixels imaging the same scene [ Fig. 2(e)]. The results from our computation are tabulated in Table 3. In addition to the estimated lifetimes and distances, we tabulate two relevant metrics.
1. First, the mean squared error or MSE, a measure of distortion, is defined as follows: let ν be the oracle estimate and fν n g N −1 n0 be N estimated values of ν. The MSE is We compute log ffiffiffiffiffiffiffiffiffi ffi MSE p forτ,d , and ∠h in Table 3, where the last term is due to fitted measurements

Research Article
which are synthesized after estimating fτ;d g.
Based on the estimates in Table 3, the expected lifetime is τ 31.3142 ns, and the estimated expected distance isd 1.0799 m, both consistent with the ground truth.

B. FD-ToF FLI: Theory and Experiments
As described in Section 2.B, in the FD-ToF mode of operation, the ToF sensor probes the scene with an AMCW of form pt 1 p 0 cosωt. Following Eq. (3), the reflected signal is The reflected signal is cross correlated [Eq. (10)] at the ToF sensor to produce measurements mt 10 jĥ0j jĥωj To estimate the phase and amplitude, we utilize the method from Section 2.B. Noting, however, that the scene transfer function is not constant, we express z explicitly as a function of the modulation frequency ω 0 : z zω 0 . The amplitude and phase estimates are (33) respectively. Note that the estimation of amplitude jĥω 0 j requires the knowledge of p 0 . However, this is not the case for phase estimation, which is prescribed by ∠ĥω 0 and is devoid of p 0 . Again, our method is blind in that we can directly estimate ∠ĥω 0 from zω 0 even if we do not know anything about p 0 .
Note in passing that multiple frequency measurements can be used to estimate multiple depths [16,20,21,31], for which ht Based on the image formation model for the SRF in Eq. (13), we will next show how multiple frequency measurements of the form fzkω 0 g K k1 can be used alternately to estimate the parameters of interest, that is, fd ; τg in the context of FLI.

Experimental Verification of FD-ToF FLI
Using the same physical setup as that in Section 3.A.2, we move the sample to d 2.5 m, and we set ω 0 ∕2π ≡ f 0 1 MHz and acquire equispaced ToF measurements fzkf 0 g K k1 ; k 1; 2; …; 40: (34) The amplitudes and phases of z for k 20, 30, and 40 are shown in Fig. 4(a). The effects of fluorescence are clearly visible in both amplitude and phase. The dielectric filter eliminates all nonfluorescent emission, leading to a noisy background signal. The fluorescing quantum dot causes an increase in the measured phase. The increased phase measurement at the sample location is due to the presence of θ τ term in Eq. (20). In context of the experimental setup described in Fig. 2(e), we mark the background pixel x b ; y b as well as the fluorescent pixel x f ; y f in Fig. 4(a). The average phase of the background pixel is noted to be ∠z x b ;y b 40f 0 4.1625 rad. This amounts to a depth of 2.4843 m, which is consistent with the experimental setup where d 2.5 m. On the other hand, the average phase at the location of the quantum dot is observed to be ∠z x f ;y f 40f 0 5.5822 rad. This relates to an erroneous depth of 3.3316 m, which is a result of multipath interference [20].
By using ToF phase measurements [Eq. (34)] for 10 different pixel locations corresponding to x f ; y f , we solve for the nonlinear least squares problem: As before in Section 3.A.2, we use the trust-region-based algorithm with the least absolute residual criterion. As a result of fitting, we estimated andτ for each pixel. The estimated distance and lifetime values together with log ffiffiffiffiffiffiffiffiffi ffi MSE p forτ,d , and ∠h for all the pixels are tabulated in Table 4. The average lifetime and distance is estimated to be fd ;τg f2.4961 m; 31.024 nsg: Both the measured phase f∠zkf 0 g K 40 k1 and the fit obtained by Eq. (35) for four pixels are plotted in Fig. 4(b).

A. ToF Sensors for Microscopy
The method generalizes to microscopy modes. A salient feature about the technique is that it is a local, per-pixel calculation, so that the lateral scale of the problem does not influence the technique. Thus, our proof-of-principle demonstration can be extended to integration with microscopy, both wide-field and pointscanning techniques, provided there is no pixel crosstalk. In fact, our current system is not aberration-corrected, and the resulting model mismatch makes reconstruction more challenging. A wellcalibrated microscopy setup should alleviate this mismatch and improve results. The important difference here is that the present approach simultaneously estimates lifetime and distance based on the same measurement. Numerical simulations in Fig. 5 indicate that our method is suitable for recovering a 4 ns lifetime. We simulate a system bandwidth of 40 MHz (half the experimental bandwidth [29]). Numerically, we compare the estimation accuracy of τ 1 4 ns and τ 2 32 ns with d 2.5 m. We do so by studying the upper bounds on the ffiffiffiffiffiffiffiffiffi ffi MSE p linked with estimation of the set of parameters of interest, fτ 1 ; dg and fτ 2 ; dg. For a fixed ToF sensor bandwidth, we vary the number of measurements [Eq. (33)] by sampling with equispaced frequencies from 0 to 40 MHz. We use four different step sizes, f k 0 (corresponding to N k samples), specified by where ε n represents additive white Gaussian noise for a SNR ranging from 0 to 60 dB. Each SNR value is averaged over 2000 realizations. In Fig. 5(a) we plot the 2000 estimated values of the lifetime parameters for each value of SNR. Intuitively, as the (a) Multifrequency measurements of τ 32 ns quantum-dot-based fluorescent sample. We show phase and amplitude images, that is, f∠z10f 0 ; ∠z20f 0 ; ∠z30f 0 ; ∠z40f 0 g (in radians, 0; 2π) and fjz10f 0 j; jz20f 0 j; jz30f 0 j; jz40f 0 jg (in decibels), respectively. The base modulation frequency for the experiment is f 0 ω 0 ∕2π 1 MHz. The phase at the background pixel is recorded to be ∠z x b ;y b 40f 0 4.1625 rad, which amounts to a depth of 2.48 m, which is consistent with the experimental setup. At the location of the fluorescent sample, we recorded a higher phase value ∠z x f ;y f 40f 0 5.5822, which is attributed to the fluorescence phenomenon. (b) Multifrequency raw phase measurements f∠z x;y kf 0 g K 40 k1 for four pixels. The measurements confirm with the theoretical hypothesis of Eq. (20), as well as the fitted phase obtained by Eq. (35). The estimated phase contribution [Eq. (20)] due to distance θd and lifetime θτ is also plotted.

Research Article
SNR increases, the estimates cluster around the oracle estimate of τ 1 4 ns and τ 2 32 ns, respectively. In Fig. 5 As the number of measurements N increases, the log ffiffiffiffiffiffiffiffiffi ffi MSE p τ drops consistently. This observation is consistent throughout all of our experiments. In fact, we note that ffiffiffiffiffiffiffiffiffi ffi MSE p τ 1 < ffiffiffiffiffiffiffiffiffi ffi MSE p τ 2 , implying that the 4 ns sample may be estimated with higher accuracy when compared to the 32 ns sample. As noted in Table 4, the operational SNR of our system is around 45 dB, implying that the lifetime parameter can be estimated with sufficient accuracy, as plotted in Fig. 5(b).
Thus, although our demonstration here operates at a lower time scale (Fig. 1) than is typical in practice, this is not a fundamental limitation of the method. The reason is that we compensate for lower time resolution by utilizing a computationally different method of inversion. Indeed, appropriate modeling and prior information lend themselves naturally to ToF sensing and offer a path toward superresolution [21,33].
Of course, calculation of the theoretical lower bounds (on lifetimes and distances) requires calculation of the Cramer-Rao bound. Though beyond the scope of the present work, this limit is ultimately dictated by noise. For instance, the Cramer-Rao lower bound for distance estimation is derived in [21] and obeys a version of the law in Eq. (36). For typical systems, we expect from Fig. 5 that the method is suitable for current application needs.

C. Experimental and Computational Precision
The current optical ToF technology allows for range estimation in millimeter precision. The variation in estimated lifetime is due to sample inhomogeneity, nonuniform lighting [34], and potential model mismatch from lens aberration.
The measurement precision here is indicated in Fig. 6, which shows cross sections of the four phase plots in Fig. 4(a). The phase contribution the first 50 pixels is due to the background. We mark the average phase value on the y axis. Using θ d f 4πf d∕c, we estimate the distance d given θ d f and f f10; 20; 30; 40g MHz. At each modulation frequency, the estimated distance is d f2.56; 2.52; 2.47; 2.48g m, whereas the actual distance is 2.5 m. The variability in the distance estimates is mainly because, across modulation frequencies, the phasefrequency relation θ d f 4πf d∕c may not be strictly linear due to distortions.

CONCLUSIONS
In conclusion, we have demonstrated a FLI alternative that is based on cost-effective ToF sensors. We simultaneously estimate the lifetime and the distance of the sample from the sensor. Unlike existing methods [8,12], our approach is calibration-free and Phase Value (in radians) requires no prior information on the experimental path length and, thus, allows for faster acquisition time. Furthermore, our technique is blind in that we do not assume the knowledge of the illumination waveform. Overall, our system shows promise for two-dimensional imaging, and can be generalized to volumes [35]. Because the technique is modular, it can be implemented with other computational imaging techniques to create a new platform for wide-field FLI. The method offers new possibilities for open questions. The case of multiple lifetime imaging [10,13] is interesting and is yet to be explored in the context of ToF sensors, both theoretically as well as experimentally. Comparison with other FLI techniques [36] and fundamental resolution limits for simultaneous estimation of lifetime and depth information will allow for better understanding of the applicability of ToF sensors for bio-imaging tasks such as tumor detection [2,3] and fluorescence tomography [4,5].