606

IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 18, NO. 2, MARCH 2014

Blind End-Member and Abundance Extraction for Multispectral Fluorescence Lifetime Imaging Microscopy Data Omar Gutierrez-Navarro, Daniel U. Campos-Delgado, Edgar R. Arce-Santana, Martin O. Mendez, and Javier A. Jo

Abstract—This paper proposes a new blind end-member and abundance extraction (BEAE) method for multispectral fluorescence lifetime imaging microscopy (m-FLIM) data. The chemometrical analysis relies on an iterative estimation of the fluorescence decay end-members and their abundances. The proposed method is based on a linear mixture model with positivity and sum-to-one restrictions on the abundances and end-members to compensate for signature variability. The synthesis procedure depends on a quadratic optimization problem, which is solved by an alternating least-squares structure over convex sets. The BEAE strategy only assumes that the number of components in the analyzed sample is known a spriori. The proposed method is first validated by using synthetic m-FLIM datasets at 15, 20, and 25 dB signal-tonoise ratios. The samples simulate the mixed response of tissue containing multiple fluorescent intensity decays. Furthermore, the results were also validated with six m-FLIM datasets from fresh postmortem human coronary atherosclerotic plaques. A quantitative evaluation of the BEAE was made against two popular techniques: minimum volume constrained nonnegative matrix factorization (MVC-NMF) and multivariate curve resolution-alternating least-squares (MCR-ALS). Our proposed method (BEAE) was able to provide more accurate estimations of the end-members: 0.32% minimum relative error and 13.82% worst-case scenario, despite different initial conditions in the iterative optimization procedure and noise effect. Meanwhile, MVC-NMF and MCR-ALS presented more variability in estimating the end-members: 0.35% and 0.34% for minimum errors and 15.31% and 13.25% in the worst-case scenarios, respectively. This tendency was also maintained for the abundances, where BEAE obtained 0.05 as the minimum absolute error and 0.12 in the worst-case scenario; MCR-ALS and MVCNMF achieved 0.04 and 0.06 for the minimum absolute errors, and 0.15 and 0.17 under the worst-case conditions, respectively. In addition, the average computation time was evaluated for the synthetic datasets, where MVC-NMF achieved the fastest time, followed by BEAE and finally MCR-ALS. Consequently, BEAE improved MVC-NMF in convergence to a local optimal solution and robustness against signal variability, and it is roughly 3.6 time faster than MCR-ALS.

Manuscript received October 3, 2012; revised June 12, 2013; August 8, 2013; accepted August 13, 2013. Date of publication August 22, 2013; date of current version March 3, 2014. The work of O. Gutierrez-Navarro was supported by a grant from the National Council of Science and Technology (CONACYT) of Mexico. This work was supported in part by CONACYT-Texas A&M University under Grant #2012-034, and in part by CONACYT under Grant #168140. O. Gutierrez-Navarro, D. U. Campos-Delgado, E. R. Arce-Santana, and M. O. Mendez are with the Facultad de Ciencias, Universidad Autonoma de San Luis Potosi, San Luis Potosi 78290, Mexico (e-mail: [email protected]; ducd@ fciencias.uaslp.mx; [email protected]; [email protected]). J. A. Jo is with the Department of Biomedical Engineering, Texas A&M University, College Station, TX 77843 USA (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/JBHI.2013.2279335

Index Terms—Autofluorescence, blind source separation, endmember extraction, fluorescence imaging, linear spectral unmixing, quadratic optimization.

NOMENCLATURE ALS AAE BEAE IC LDL MCR-ALS

Alternating least-squares. Abundance average error. Blind end-member and abundance estimation. Initial condition. Low-density lipoproteins. Multivariate curve resolution-alternating least-squares. m-FLIM Multispectral fluorescence lifetime imaging. NMF Nonnegative matrix factorization. MVC-NMF Minimum volume constrained nonnegative matrix factorization. PRE Profiles relative error. SE Synthetic end-member. SNR Signal-to-noise ratio. I. INTRODUCTION EVERAL research groups around the world study options for an effective method to perform minimal invasive characterization of living tissue [1]. Fluorescent dyes are a common tool employed to label molecules within living systems, since this technique allows their location and therefore a description of their molecular environment, which could be useful for pathology detection and treatment monitoring in a noninvasive way at low cost [2]. The estimation of the fluorescent components and their fractional contributions is known in the literature as the linear unmixing problem [3]. Previous approaches look to achieve this objective by using synthetic dyers [4], [5], which rely on invasive techniques for their application. m-FLIM is a recently developed technique [6] which measures the simultaneous time-resolved fluorescence response at different wavelength bands of the molecular content of a tissue sample. Fig. 1 shows an endogenous m-FLIM measurement performed on a postmortem human coronary artery sample, where the fluorescence emission was measured at three distinctive wavelength bands (390 ± 20, 452 ± 22.5, and 550 ± 20 nm). m-FLIM data take advantage of the autofluorescence decay profiles to reduce the need for synthetic dyers, although its main limitation is the similarity among the autofluorescence profiles or end-members from different molecules that build up the

S

2168-2194 © 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications standards/publications/rights/index.html for more information.

GUTIERREZ-NAVARRO et al.: BLIND END-MEMBER AND ABUNDANCE EXTRACTION

607

Fig. 1. m-FLIM data arranged in a 3-D stack, which can be seen as an array of images. In the experiments, the data have a dimension of 60 × 60 × 510 samples. Each measurement in the spatial plane (x, y) contains the concatenated fluorescence decay measurements from three wavelength bands, namely 390 ± 20, 452 ± 22.5, and 550 ± 20 nm.

Fig. 2. Plot of three endogenous fluorescence end-members corresponding to collagen, elastin, and LDL. These end-members contain the concatenated fluorescence decays measured at 390 ± 20, 452 ± 22.5, and 550 ± 20 nm. Although these end-members are generated by different molecular components, they present high overlapping among them.

tissue sample. Consequently, there is high overlapping between the end-members in the recorded measurements. This effect is illustrated in Fig. 2. Hence, advanced blind source separation techniques are needed to effectively characterize the m-FLIM samples [7], [8]. The linear unmixing problem is found in a large number of applications such as remote sensing [9], [10], where techniques such as nonnegative matrix factorization have been extensively used. Most of the existing techniques rely on an optimization characterization to solve the unmixing problem [9]. In fact, sev-

eral techniques are based on the fact that due to the linear mixture model, end-members represent the vertices of an N -dimensional simplex [11], where N is the number of components in the mixture. Different strategies are employed to estimate these vertices; for example, one approach searches for a minimumvolume simplex in order to fit the measured data [12], [13]. However, the resulting optimization strategy is nonlinear and its convergence depends largely on the ICs. Another strategy relies on the Bayesian estimation for joint end-member extraction and abundance estimation in hyperspectral imagery, where a Gibbs sampler is proposed to reconstruct the posterior distributions [14]. However, the complexity of the resulting solution could be a limiting factor for large data samples. Meanwhile, in [7], a quadratic cost function is proposed which avoids the estimation of the simplex volume for end-member extraction in hyperspectral imagery, and has the advantage of a global solution. In this context, this paper considers the quadratic formulation in [7], [15], and [16], and extends them by adding new restrictions particularly of m-FLIM data for BEAE. The mathematical derivations of the new methodology are included in the paper, and the convergence of the BEAE algorithm is analytically demonstrated based on ALS. Synthetic data are used to systematically evaluate the performance and robustness of the BEAE, in comparison with two of the most common accepted techniques in the literature: MVC-NMF and MCR-ALS. Finally, six m-FLIM datasets from fresh postmortem human coronary atherosclerotic plaques are also employed to validate our proposal. The rest of the paper is organized as follows. The linear mixture model for m-FLIM data and the notation used throughout the paper are presented in Section II. Section III describes the proposed methodology based on quadratic optimization and ALS. The synthetic evaluation for end-members and abun-

608

IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 18, NO. 2, MARCH 2014

dances estimation are described in Section IV. In Section V, the results obtained with ex vivo samples are detailed, and finally, the discussion of the results and conclusions are presented in Section VI. The notation used in this paper is described next. Scalars are denoted by italic letters, and vectors and matrices by lower and upper case bold letters, respectively. R represents the real numbers, and card(Ω) the cardinality of the set Ω. For a real vector x, the transpose√operation is denoted by x , the Euclidean norm by x2 = x x, and x  0 represents that each component in the vector is positive or zero. For a matrix X, Xi,jrepresents the element in the i-row and j-column,Tr(X) = i Xi,i denotes the trace operation, and XF = Tr (X X) the Frobenius norm. X  0 represents that each component in the matrix is nonnegative. An N-dimensional vector filled with ones is represented by 1N , and IN denotes the identity matrix of order N . II. LINEAR MIXTURE MODEL FOR M-FLIM DATA According to the model in [3], every fluorescence decay mixture yk ,λ (t) is a linear combination of the fluorescence decay profiles or end-members pn ,λ (t), and their fractional contributions or abundances αk ,n . From now on, we will refer to the fluorescence decay profiles simply as profiles or end-members. At a certain spatial position k and wavelength λ, over t = 1, . . . , T time samples, the data mixture is modeled as yk ,λ (t) =

N 

αk ,n pn ,λ (t)

∀k = 1, . . . , K, & λ = 1, . . . , W

n =1

(1) where the abundance coefficients αk ,n are normalized at each  spatial position k, i.e., N n =1 αk ,n = 1, K is the number of spatial locations in the sample, N is the number of end-members, and W is the number of wavelength bands. The goal of blind unmixing is to estimate the end-members pn ,λ , and their abundances αk ,n , which compose the mixture yk ,λ sampled at every spatial position k of the m-FLIM data. The mixture model in (1) can be written in a vector notation for the m-FLIM data over the sampled time-window as yk ,λ =

N 

αk ,n pn ,λ

(2)

n =1

where yk ,λ = [yk ,λ (1) · · · yk ,λ (T )] and pn ,λ = [pn ,λ (1) · · · pn ,λ (T )] are the measurement vector at k space sample and the end-member vector for nth component and λ wavelength, respectively. To further simplify the notation, the information at every position k is concatenated for all wavelength bands λ = 1, . . . , W as a single vector for the nth end-member:   T ∈ RD , ∀n = 1, . . . , N pn = pTn ,1 , . . . , p n ,λ , . . . , pn ,W (3) where D = W · T is the total number of time measurements gathering all wavelength bands, W is the number of wavelength bands, and T is the number of time samples. The mixture samples are formatted in a similar fashion   ∈ RD ∀k = 1, . . . , K. (4) yk = yk,1 , . . . , yk,W

Hence, the linear model in (1) can be written in a matrix notation by collecting the whole time samples and wavelengths as yk = [p1 , · · · , pN ] [αk ,1 , · · · , αk ,N ]  

 

P∈RD ×N

∀k = 1, . . . , K. (5)

αk ∈RN

Finally, by incorporating all the spatial samples, the mixture model can be described by three matrices: Y = [y1 , . . . , yK ] = [p1 , . . . , pN ] [α1 , . . . , αK ]  

 

(6)

A

P

where A ∈ RN ×K is a matrix containing the abundances of every end-member at each available position of the sample data, and Y ∈ RD ×K is the overall m-FLIM measurement matrix. To avoid redundancy in the linear mixture model in (2), the end-members (p1 , . . . , pN ) are assumed linearly independent and D  N , and as a result, P is a full-column rank matrix. Meanwhile, the abundance matrix A is assumed full-row rank, since it is assumed K  N . In general, when working with fluorescence emissions, it is difficult to measure the exact same intensity from a sample, even if it contains the same components and it is under controlled conditions, and specially from different imaging sessions. This effect is due to many artifacts that may arise in fluorescence imaging devices [17], [18]; this phenomenon is known as signature variability in the hyperspectral imaging literature [19]. These variations in the fluorescence intensity measurements will result in differences between the decay profiles estimated in different samples even though they contain the same fluorophores and remain with the same concentrations. Another effect due to these fluctuations is that the mixtures might not be contained within the simplex described by the original components, which violates constraint (8) [20], [21]. To minimize the effects of signature variability, it is a common practice to impose a normalization constraint on the measurement data at each spatial location k, so that it is satisfied yk 1D = 1∀ k = 1, . . . , K. As a consequence, the end-members are also with reT  normalized p spect to time and wavelength bands, i.e., W λ=1 t=1 n ,λ (t) = 1 ∀n = 1, . . . , N . Furthermore, every intensity value in the endmembers matrix P has to be nonnegative. Therefore, considering all the restrictions and measurements available, the complete linear unmixing problem for m-FLIM data can be written as an approximation problem [9] 1 min Y − PA2F P,A 2

(7)

such that α k 1N =

N 

αk ,n = 1

∀k

(8)

n =1

αk ,n ≥ 0 p n 1D =

∀k, n

T W  

pn ,λ (t) = 1

(9) ∀n

(10)

λ=1 t=1

pn ,λ (t) ≥ 0

∀n, λ, t.

(11)

GUTIERREZ-NAVARRO et al.: BLIND END-MEMBER AND ABUNDANCE EXTRACTION

Consequently, the cost function in (7) evaluates the similarity between the measurements and predicted linear mixtures. In fact, the constraints in the vectors pn and αk can be expressed as conditions in the matrices A and P, which can define the following feasible sets:  (12) ΩA = A ∈ RN ×K : A 1N = 1K and A  0  ΩP = P ∈ RD ×N : P 1D = 1N and P  0 . (13) From these definitions, the feasible sets ΩA and Ωp are convex (see Appendix A), and thus, the quadratic approximation problem in (7) is solved over convex feasible sets. This property will be crucial in order to guarantee a global optimum solution at each step of our iterative search regardless of the ICs, as will be shown in the following section.

609

regularization term derived in Appendix B is substituted:

 

1 J = Y − PA2F + χ P 1D − 1N + ρTr POP 2 (15) where χ ∈ RN is a vector of Lagrange multipliers and O = N IN − 1N 1 N is a symmetric matrix (see Appendix C). The optimization of the objective function J in (15) implies the following stationary conditions:

 ∂J = −YA + P AA + ρO + 1D χ = 0 ∂P ∂J = P 1D − 1N = 0. ∂χ

III. BEAE FOR M-FLIM DATA

A. Quadratically Constrained Blind End-Member Extraction Given N components in the sample, and their fractional contributions A, our methodology focuses on the following quadratic optimization problem to obtain the profiles matrix P: min

P∈Ω P

N −1  N  1 Y − PA2F + ρ pi − pj 22 2 i=1 j =i+1  

(17)

Furthermore, it is satisfied ∂2 J = AA + ρO > 0, ∂P2

In this study, our previous result on the abundance estimation [15], [16] is employed in conjunction with a new strategy to estimate the end-members. The methodology assumes that the number of components in the sample is known, and consequently, it is needed to estimate the end-members and their respective fractional contributions or abundances. Moreover, this strategy relies on ALS [22] in order to iteratively compute the blind end-member extraction, and next the abundances estimation. This new scheme is based on the Berman et al. [7] algorithm, but extended to include the sum-to-one and positivity restrictions in the resulting end-members.

(16)

(18)

since AA is positive definite, and O is symmetric and diagonally dominant with positive entries in the main diagonal, so O is a positive semidefinite matrix [25], which guarantees a minimum for the solution of the stationary conditions in (16)–(17). The optimal closed-form solution is then given by  

−1 1 1  + 1D 1 P = ID − 1D 1D YA AA + ρO D D D (19) where the mathematical derivations are detailed in Appendix C. In fact, the solution in (19) could produce negative elements, since only the sum-to-one condition P 1D = 1N is enforced by the Lagrange multiplier method. As a result, an orthogonal projection is then proposed to the feasible set ΩP for any negative elements in (19) [24], i.e., any negative element in (19) is set to zero and the rest of the elements each per column of P are normalized. B. Abundance Estimation

(14)

L (P,A)

 −1 n where the quadratic regularization term ni=1 j =i+1 pi − 2 pj 2 proposed in [7] is included in order to avoid an ill-posed optimization problem [9], and ρ > 0 is a weighting factor to control the influence of the regularization term. This term is related to the sum of the distances between the vertices of the simplex generated by the end-members (p1 , . . . , pN ). These distances should be minimized in order to search for the simplex with the minimum volume containing the measured data [7]. Contrary to other approaches in [5], [12], and [23], the quadratic regularization term simplifies the overall optimization process in order to reach for a closed-form solution. In fact, from the synthesis formulation in (14), a quadratic optimization problem is proposed over a convex feasible set ΩP , and as a result, an optimum solution can always be guaranteed. Next, the cost function L(P, A) in (14) is extended by the equality restrictions P 1D = 1N (end-members normalization) through Lagrange multipliers [24], and also the equivalent representation of the

Once there is an estimation on the end-members P, the second stage is to calculate the abundances matrix A according to the following constrained quadratic optimization problem [16]: min

A∈Ω A

1 Y − PA2F . 2

(20)

The previous overdetermined optimization is solved over the convex set ΩA , which guarantees a global solution. Since the abundances are just dependent on the spatial location, the optimization in (20) can be performed at each spatial measurement k: 1 min yk − Pαk 2 = αk 2  2 W N   1    min αk ,n pn ,λ  (21) yk ,λ − α k , 1 ,...,α k , N 2   λ=1

n =1

2

such that α k 1N = 1 αk  0.

(22)

610

IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 18, NO. 2, MARCH 2014

Fig. 3. Difference of the cost function L i for successive iterations. The decreasing pattern of | L(P i + 1 , Ai + 1 ) − L(P i , Ai ) | illustrates the local convergence of the BEAE algorithm.

atively to construct a solution sequence (Pi , Ai ) such that

 Ai+1 = arg min L Pi , A A∈Ω A

i+1

P



= arg min L P, Ai+1 P∈Ω P

(27)

where it is then satisfied



 L Pi , Ai+1 ≤ L Pi , Ai



 L Pi+1 , Ai+1 ≤ L Pi , Ai+1

Therefore, we consider our previous closed-form solution in [16] and [15] that is described in Algorithm 1. This algorithm requires prior knowledge on the number of components N and the endmembers matrix P. A more detailed information regarding the matrices employed in the algorithm can be found in Appendix D. In fact, a quadratic regularization term is not needed in (21), as in (14), since this optimization scheme is not ill-posed for linearly independent end-members (p1 , . . . , pN ) [15], [16], i.e., there is a unique solution at each spatial location [24].

due to the quadratic dependence of L on the end-members and abundance matrices (P, A). Therefore, by defining 

(28) Li+1 = L Pi+1 , Ai+1 , a decrement after each complete iteration is guaranteed: Li+1 ≤ Li

∀i ≥ 0.

(29)

Fig. 3 shows an example of this decreasing pattern for one of the synthetic experiments described in the following section. Furthermore, since the cost function L has a continuous dependence on the matrices (P, A), the local convergence of (27) is also guaranteed. The overall iterative methodology is defined in Algorithm 2.

C. BEAE for m-FLIM Data First, it is important to point out that the cost function L(P, A) in (14) depends on the product of two unknowns: end-members P and abundances matrix A. Consequently, in order to reach a solution to the joint minimization problem in (14), an ALS formulation is adopted [22]. In this way, the solution to the quadratically constrained blind end-member extraction in (19), and the abundance estimation of Algorithm 1, are applied iter-

IV. SYNTHETIC EVALUATION In this section, we generate synthetic m-FLIM measurements [15], [16] in order to provide a controlled scenario for performance tests. The proposed method, called BEAE, was evaluated and compared against two well-known iterative techniques: MVC-NMF [23] and MCR-ALS [26], which can deal with the nonnegativity constraint in (9). These methods are also

GUTIERREZ-NAVARRO et al.: BLIND END-MEMBER AND ABUNDANCE EXTRACTION

611

Fig. 4. Five end-members employed in the synthetic simulations for performance evaluation. The similarity among the SEs is high to simulate real endogenous components. The plot uses a semilog scale in the y-axis to enhance the differences among the amplitudes of the SEs.

capable of unmixing the m-FLIM data into multiple components, and require an initial estimation of the number of components, as well as initial end-member candidates, which will be referred as IC. Two ICs were employed: the first one was chosen by spatial positions corresponding to pure samples of each end-member and the second one was set to initial candidates corresponding to mixtures of the pure samples that highlights a more difficult scenario. In order to provide a fair comparison, the sum-to-one constraint in (11) was imposed by normalizing the end-members at the end of each iteration in MVC-NMF and MCR-ALS. Meanwhile, the sum-to-one constraint for the abundances in (8) was also forced by a proper projection and scaling, since these methods do not satisfy it strictly. Synthetic data were first generated in order to evaluate the performance of the blind end-member extraction algorithms. Two scenarios were conducted by producing synthetic data with three and five SEs. However, since both evaluations provided comparable quantitative results, just the test with the larger number of components is illustrated next. A total of five SEs were constructed by using the stretched-exponential function from [15] and [27] to simulate the intensity decay from endogenous fluorophores. The resulting end-members are depicted in Fig. 4. They were employed as ideal end-members to generate the synthetic m-FLIM data. These end-members, referred as true profiles from now on, present high overlapping, and consequently, the proposed unmixing problem presents important numerical challenges. An abundance map was generated to simulate a highly mixed scenario without the presence of pure samples from a single component. Additive Gaussian noise was added to characterize the sensitivity of the m-FLIM data under the linear unmixing model as Y = PA + Υ

(30)

where Υ is a noise matrix of independently distributed Gaussian samples with zero mean. The variance of the noise elements was set to satisfy 15, 20, and 25 dB SNRs. The evaluated algorithms were initialized by using mixture samples randomly selected from the evaluated synthetic m-FLIM datasets. Two different ICs for the end-members matrix P0 were selected to evaluate the sensitivity of the solutions with respect to the initial point in the iterative schemes. For each IC and SNR, 30 m-FLIM datasets were considered for each algorithm. In the evaluation of BEAE, the weight on the regularization parameter ρ in (14) was set to 0.05, and the stoppage tolerance  to 0.001. A. End-Member Extraction The estimated end-members were compared against the true fluorescence profiles from Fig. 4. The relative error was employed to asses quantitatively the quality of the extracted endmembers from the synthetic m-FLIM data. The PRE for a given component n is defined as PREn = 100 ×

 n 2 pn − p pn 2

∀n = 1, . . . , N,

(31)

 n denotes the estimated end-member and pn the true where p vector for the nth component. The mean and standard variation of the PRE were estimated for each experiment comprising 30 simulations, and the results for each SE are shown in Table I. In this table, the value in boldface per row indicates the best score, or values that are within 5% tolerance of this best result. The mean PRE for BEAE is the lowest in the majority of experiments (70% of cases) under a large noise or uncertainty environment (15 dB), and also the PRE presented the lowest standard deviation (lowest variability) under this same condition, as noted in Table I. In the meantime, under low noise or uncertainty conditions (25 dB), all three algorithms reported similar mean PRE performance, but in almost all cases BEAE reported less variability (standard deviation). The results from

612

IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 18, NO. 2, MARCH 2014

TABLE I MEAN ± STANDARD DEVIATION OF THE PRE BETWEEN THE TRUE AND ESTIMATED END-MEMBERS (THE LOWEST MEAN VALUES PER SCENARIO ARE HIGHLIGHTED WITH BOLDFACE)

Table I also show that BEAE obtained for its best estimation a 0.32% relative error, meanwhile 13.82% error for the worst-case condition, despite the ICs in the iterative optimization procedure plus noise effects. Meanwhile, MVC-NMF and MCR-ALS reported 0.25% and 0.24% minimum mean PRE, and 15.31% and 13.26% in their worst-case performances, respectively. An important feature of MVC-NMF and MCR-ALS is that they consistently presented more variability in estimating the endmembers compared to BEAE, as can be seen from the standard deviations of the PRE. B. Abundance Estimation The abundances were evaluated in a similar fashion to the endmembers. The estimated fractional contributions were compared by using the abundance average error (AAE) computed over the whole sample set as   K 1  ( αk ,n − αk ,n )2 , ∀ n = 1, . . . , N AAEn =  K k =1

(32) where α k ,n are the estimated abundance and αk ,n the true value for k spatial location and n end-member. For illustration purposes, the abundance maps from one experiment with the SNR equal to 15 dB are shown in Fig. 5. Table II shows the mean and standard deviation values of the AAE measured for each algorithm under different ICs and noise levels. Once again, the values in the boldface represent the lowest values in each testing scenario with a range of tolerance of 5%. One can observe from the table that MCR-ALS and BEAE

TABLE II MEAN ± STANDARD DEVIATION OF THE AAE BETWEEN THE TRUE AND ESTIMATED ABUNDANCE MAPS (THE LOWEST MEAN VALUES PER SCENARIO ARE HIGHLIGHTED WITH BOLDFACE)

presented the lowest mean values at different ICs and SNR, but the former achieved a higher variance in the AAE. BEAE reported more accurate results in the presence of large noise or uncertainty (15 dB). Moreover, BEAE obtained a mean AAE between 0.0507 and 0.1163 for all conditions with very low variability with respect to the other two methods. In the meantime, MVC-NMF and MCR-ALS achieved 0.0648 and 0.0415 for its minimum mean AAE, and 0.1679 and 0.1484 under the worstcase conditions, respectively. These results illustrate that BEAE presents more accurate estimation independent of the SNR and ICs, resulting in a more robust approach. Meanwhile, MVCNMF achieves the worst-error performance with high variability for all scenarios. V. EXPERIMENTAL EVALUATION The blind end-member extraction algorithms (BEAE, MVCNMF, and MCR-ALS) were also tested by using m-FLIM datasets from fresh postmortem human coronary atherosclerotic plaques [6], [28]. Nonetheless, there is no quantitative description available that could be used as ground truth for a clean comparison, so only qualitative evaluations are provided in this section. In this way, using the fluorescence intensity decays for collagen, elastin, and LDL as ICs for the end-members, the solution provided by the three algorithms is consistent, as can be seen in Fig. 6(a), where one sample containing mostly elastin and LDL was decomposed by using BEAE, MVC-NMF, and MCR-ALS. On the other hand, when a mixture of the ideal endmembers were employed as an IC, only BEAE and MCR-ALS

GUTIERREZ-NAVARRO et al.: BLIND END-MEMBER AND ABUNDANCE EXTRACTION

613

Fig. 5. Simulated abundance map (noise free), and the results obtained using BEAE, MVC-NMF, and MCR-ALS algorithms. These results were obtained from a single experiment by using IC #2 and SNR = 15 dB.

Fig. 7. Estimated end-member profiles resulted from applying the BEAE algorithm to six m-FLIM datasets from fresh postmortem human coronary atherosclerotic plaques.

Fig. 6. Estimated abundance maps from an ex vivo sample containing mostly elastin and LDL. BEAE, MVC-NMF, and MCR-ALS were employed to unmix the sample by using different ICs, where MVC-NMF was the most sensible algorithm to changes in this parameter. (a) IC: ideal end-members. (b) IC: mixed end-members.

provided similar solutions to the one obtained previously, but not MVC-NMF, as can be seen in Fig. 6(b). This phenomenon reflects the sensitivity of MVC-NMF to the ICs. We omit more ex vivo results with MVC-NMF and MCR-ALS, since a quantitative evaluation cannot be achieved. BEAE was further validated by using a total of six m-FLIM ex vivo datasets [6], [28]. The m-FLIM data from the six samples were pooled together to estimate both the end-member and component abundances. The estimated end-members are shown in Fig. 7. The profiles of the first component (Fig. 7, top panel) showed normalized intensities of 41%, 32%, and 27% in the first, second, and third spectral channels, respectively, plus

614

IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 18, NO. 2, MARCH 2014

Fig. 8. Representative histopathology slides (left panels) and estimated abundance maps (right panels) from three coronary samples showing relatively uniform histopathology evaluation.

Fig. 9. Representative histopathology slides (left panels) and estimated abundance maps (right panels) from three coronary samples showing mixed histopathology evaluation.

average lifetime values between ≈ 7 and 11 ns, resembling the fluorescence emission of LDL. The profile of the second component (see Fig. 7, middle panel) showed normalized intensities of 45%, 37%, and 18% in the first, second, and third spectral channels, respectively. The estimated average lifetime values were between ≈ 7 and 8 ns, which are similar to the fluorescence emission of collagen. Finally, the third component (see Fig. 7, bottom panel) showed normalized intensities of 45%, 37%, and 18% in the first, second, and third spectral channels, respectively, and average lifetime values between ≈ 6 and 7 ns, resembling the fluorescence emission of elastin. Thus, the first, second, and third estimated end-member profiles and the corresponding abundance maps were associated with LDL (lipids), collagen, and elastin, respectively. Three of the six coronary samples presented a relatively uniform histopathology, as illustrated in Fig. 8. The first sample (see Fig. 8, top panel) corresponded to an almost normal coronary artery rich in elastin, shown as dark fibers in the histopathology evaluation, and the associated abundance maps showed dominant contribution from elastin (≥ 60%), with lesser contribution from collagen (0–30%) and virtually no contribution from LDL. The second sample (see Fig. 8, middle panel) cor-

responded to a thin-cap fibroatheroma with a large lipid core, shown as empty spaces in the histopathology evaluation, where the associated abundance maps presented dominant contribution from LDL (≥ 80%), with little to none contributions from collagen and elastin. The third sample (see Fig. 8, bottom panel) corresponded to a fibrotic plaque composed of mostly collagen, shown as light-blue fibers in the histopathology evaluation, and the associated abundance maps showed dominant contribution from collagen (> 90%), with virtually no contributions from neither elastin nor LDL. On the other hand, the remaining three coronary samples showed mixed histopathology samples, as shown in Fig. 9. The first sample (see Fig. 9, top panel) corresponded to a plaque with areas of significant content of elastin, collagen, and/or LDL, as shown in the histopathology slide (corresponding to the middle cross-section of the imaged plaque shown as a dotted black line in the abundance maps). The abundance maps also indicated areas dominated by LDL, collagen or elastin, reflecting a complex underlying histopathology. The second sample (see Fig. 9, middle panel) corresponded to a plaque with a thick region rich in lipids (left side on histopathology), and a thin region rich in elastin (right side on histopathology). The abundance maps also

GUTIERREZ-NAVARRO et al.: BLIND END-MEMBER AND ABUNDANCE EXTRACTION

clearly indicated two distinct regions, one dominated by LDL (top-left) and the other dominated by elastin (bottom-right), with little to none contribution from collagen. The third sample (see Fig. 9, bottom panel) corresponded to a plaque with a thin region rich in elastin (left side on histopathology) and a thick region rich in collagen (right side on histopathology). The abundance maps evidently indicated two distinct regions: one dominated by elastin (top) and the other dominated by collagen (bottom), with little to none contribution from LDL. VI. DISCUSSION AND CONCLUSION This study proposed a novel method to automatically detect fluorescence end-members and their relative abundances in a biological sample based on the analysis of m-FLIM data. A comparison of the proposed method (BEAE) with two popular methods from the literature (MVC-NMF and MCR-ALS) used for blind spectral unmixing is also presented. Our main observations were the following. 1) The experimental results support the hypothesis that fluorescence measurements could be used to quantitatively characterize a tissue sample, as a function of its endmembers and its respective abundances. 2) The BEAE methodology can be successfully employed to perform end-member extraction and abundance quantification in m-FLIM data, regardless of the number of components, and without the influence of the ICs in the iterative search. The synthetic simulations showed that BEAE is both robust to noisy or uncertain data and is tolerant to high overlapping among the fluorescence decay profiles. By comparing to other iterative schemes as MVC-NMF [23] and MCR-ALS [26], one major advantage of the BEAE algorithm is that, due to its quadratic formulation over convex feasible sets, it always converges to a minimum in the feasible region under the constraints in (8)–(10). The comparison methods evaluated in this study, MVC-NMF and MCR-ALS, are usually called blind source separation methods, although they do require some prior information in order to provide a good solution. They need to know the number of components to be solved, as well as some initial solution within the feasible region subject to the constraints (8)– (11), and close to the optimum solution for better performance. The synthetic experiments were carried out by using ICs inside the same m-FLIM datasets. Our evaluations showed that the performance of MVC-NMF could be seriously affected by the selected IC or end-members; this phenomenon is specially evident in the experimental evaluation in Fig. 6(b). As can be seen in Tables I and II, BEAE scores the smallest error in endmember and abundance estimation, when there is large noise or uncertainty in the data. MCR-ALS also obtained good performance specially when the SNR is high. However, this method has higher computational cost among the algorithms tested. The synthetic m-FLIM datasets have a dimension of 60 × 60 × 510 data, and the blind end-member extraction algorithms were implemented in MATLAB, under an operating system Ubuntu 12.04 AMD64. The computing unit has an Intel(R) Core(TM) i5-2430 CPU at 2.40 GHz with 8 GB of RAM. The average

615

run-times for MVC-NMF and BEAE were 4.9 and 8.5 s, respectively, while the average time for MCR-ALS was 30.5 s. In this way, MVC-NMF had the best processing time but the worst performance (see Fig. 6). Hence, BEAE establishes a good compromise between accuracy and computational time. In our analysis, further experimentation is needed to estimate the optimal value for the regularization term ρ in (14), and the tolerance  employed in the convergence criteria for BEAE, since both values were tuned by trial and error to have the best performance. In addition, once the end-members are estimated, the physical identification of these components from experimental databases is an another open problem. Finally, by using postmortem human coronary atherosclerotic plaques as validation elements, we have demonstrated that the application of BEAE to the experimental m-FLIM data from coronary plaques allowed estimating and identifying the three main sources of autofluorescence, namely elastin, collagen, and LDL. Furthermore, the associated abundance maps allowed the characterization of the biochemical composition of the imaged atherosclerotic plaques. These results thus demonstrate the potential application of BEAE to autofluorescence m-FLIM data for molecular imaging of biological tissue. APPENDIX A The proof of convexity of the sets ΩA and ΩP is analogous, so the discussion will be focused just on ΩA . First, take any two elements B, C ∈ ΩA and ζ ∈ [0, 1]; then, a new matrix is constructed: D = ζB + (1 − ζ)C ∈ RN ×K .

(33)

Since 0 ≤ ζ ≤ 1 and B, C  0, D  0 is satisfied. Furthermore, it is deduced that   D 1N = ζB  + (1 − ζ)C  1N = ζ1N + (1 − ζ)1N = 1N . (34) Therefore, it is concluded that D ∈ ΩA for any ζ ∈ [0, 1], and consequently, ΩA is a convex set.  APPENDIX B In this appendix, the regularization term used in the cost function (14) is simplified in order to derive an equivalent representation in terms of matrix P. From (14), and following the steps in [7], it is concluded that N −1 

N 

pi − pj 22

i=1 j =i+1

=

N  N  1 i=1 j =1

=

2

pi − pj 22

N  N  1 i=1 j =1

=N

N  i=1

2

pi 22 + pj 22 − 2p i pj

pi 22 −

N  N  i=1 j =1

p i pj



616

IEEE JOURNAL OF BIOMEDICAL AND HEALTH INFORMATICS, VOL. 18, NO. 2, MARCH 2014

=N

N 

 pi 22 −

i=1

N 

⎞  ⎛ N  ⎝ pi pj ⎠

i=1

APPENDIX D With respect to the notation described in Section II, the matrices involved in Algorithm 1 are detailed in this section: ⎡ ⎤ p1,λ , p1,λ  . . . p1,λ , pN ,λ  W ⎢ ⎥ .. .. .. A= ⎣ ⎦ (48) . . . λ=1 pN ,λ , p1,λ  . . . pN ,λ , pN ,λ 

j =1

 = N Tr PP − (P1N ) (P1N ) 



 = N Tr PP − Tr P1N 1 NP ⎞ ⎛  ⎟ ⎜  = Tr ⎝P N IN − 1N 1 N P ⎠ 

 O ∈RN ×N

(35) 

APPENDIX C



P 1D − 1N = 0.

(AA + ρO) (AY −

χ1 D)

(50)

· 1D = 1N

 AY 1D − χ1 D 1D = (AA + ρO) · 1N 1  AY 1D − (AA + ρO) · 1N ∴χ= D 1  ⇒ χ = {1 YA − 1 N (AA + ρO)} D D

λ=1

(41)

where AA + ρO is a positive definite matrix, and consequently nonsingular, and substituting into (40), 

yk ,λ , yk ,λ  , c = [1 . . . 1] 1×N , d = 1.

where μ > 0 is the Lagrange multiplier associated with the abundances sum-to-one constraint in (8). For further details, the reader is referred to [15] and [16]. 

(39)

From equation (39), matrix P can be obtained: P = (YA − 1D χ )(AA + ρO)−1 ,

W 

(40)

  1 Y − PA2F + χ P 1D − 1N + ρTr POP . 2 (36) The stationary conditions for the optimum are ∂J =0 (37) ∂P ∂J =0 (38) ∂χ and the following equations are derived: J=

  − YA + 1 D χ + P(AA + ρO) = 0

ek =

The operator ·, · denotes the inner product between vectors of the same dimensions. It is assumed that matrix A is nonsingular, i.e., the set of end-members (p1,λ , . . . , pN ,λ ) is linearly independent for each wavelength λ. In case there are negative elements in xk , the solution provided by (23), an active constraints method is employed to incorporate new equality constraints. This leads to the next system of linear equations: ⎤ ⎡ ⎤ ⎡ ⎡ ⎤ xk 2A c H . . . H −bk 1 L  0 0 ··· 0 ⎥ ⎢ μ ⎥ ⎢ d ⎥ ⎢c ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ 0 ··· 0 ⎥ ⎢ δ1 ⎥ = ⎢ 0 ⎥ ⎢ H1 0 ⎥ ⎥ ⎢ ⎢ ⎢ . ⎥ . . . . . .. ⎣ .. ⎦ (51) ⎣ .. .. .. .. ⎦ ⎣ .. ⎦ . 0 0 ··· 0 δL HL 0  

 

 

β Γ Δ

The proposed quadratic cost function in (15) is given by

−1

(49)

λ=1



= Tr POP .



W 

[yk ,λ , p1,λ  . . . yk ,λ , pN ,λ ]

bk = −2

(42) (43) (44) (45)

where 1 D 1D = D. Substituting the Lagrange multipliers vector χ in (41), we can reach the optimal solution in (19):   1 1    1 P = YA − 1D 1 YA + 1 (AA + ρO) D D D D D × (AA + ρO)−1 (46)   1 1   −1 + 1D 1 = ID − 1D 1 D YA (AA + ρO) D. D D (47) 

ACKNOWLEDGMENT The authors would like to thank Dr. H. Qi for providing the code for the MVC-NMF algorithm. They would also like to thank the anonymous reviewers who provided helpful comments to improve the quality of this study.

REFERENCES [1] M. O’Donnell, E. R. McVeigh, H. W. Strauss, A. Tanaka, B. E. Bouma, G. J. Tearney, M. A. Guttman, and E. V. Garcia, “Multimodality cardiovascular molecular imaging technology,” J. Nucl. Med., vol. 51, pp. 38S– 50S, May 2010. [2] N. Anthony, K. Berland, and P. Guo, “Principles of fluorescence for quantitative fluorescence microscopy,” in FLIM Microscopy in Biology and Medicine. London, U.K.: Chapman and Hall, Jul. 2009, pp. 35–63. [3] J. R. Lakowicz, Principles of Fluorescence Spectroscopy. New York, NY, USA: Springer-Verlag, 2006. [4] A.-S. Montcuquet, L. Herve, F. Navarro, J.-M. Dinten, and J. Mars, “In vivo fluorescence spectra unmixing and autofluorescence removal by sparse nonnegative matrix factorization,” IEEE Trans. Biomed. Eng., vol. 58, no. 9, pp. 2554–2565, Sep. 2011. [5] H. Xu and B. W. Rice., “In-vivo fluorescence imaging with a multivariate curve resolution spectral unmixing technique,” J. Biomed. Opt., vol. 14, no. 6, pp. 064011/1–064011/9, 2009.

GUTIERREZ-NAVARRO et al.: BLIND END-MEMBER AND ABUNDANCE EXTRACTION

[6] J. Park, P. Pande, S. Shrestha, F. Clubb, B. E. Applegate, and J. A. Jo, “Biochemical characterization of atherosclerotic plaques by endogenous multispectral fluorescence lifetime imaging microscopy,” Atherosclerosis, vol. 220, no. 2, pp. 394–401, 2012. [7] M. Berman, H. Kiiveri, R. Lagerstrom, A. Ernst, R. Dunne, and J. Huntington, “ICE: A statistical approach to identifying endmembers in hyperspectral images,” IEEE Trans. Geosci. Remote Sens., vol. 42, no. 10, pp. 2085–2095, Oct. 2004. [8] F. J. Theis, R. Neher, and A. Zeug, “Blind decomposition of spectral imaging microscopy: A study on artificial and real test data,” in Proc. 8th Int. Conf. Ind. Compon. Anal. Signal Separation, 2009, pp. 548–556. [9] J. M. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P. Gader, and J. Chanussot, “Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches,” IEEE J. Sel. Top. Appl. Earth Observat. Remote Sens., vol. 5, no. 2, pp. 354–379, Apr. 2012. [10] A. Plaza, G. Martin, J. Plaza, M. Zortea, and S. Sanchez, “Recent developments in endmember extraction and spectral unmixing,” in Optical Remote Sensing., (Series Augmented Vision and Reality), vol. 3. Berlin, Germany: Springer, 2011, pp. 235–267. [11] M. Craig, “Minimum volume transforms for remotely sensed data,” IEEE Trans. Geosci. Remote Sens., vol. 32, no. 3, pp. 542–552, May 1994. [12] M. B. Lopes, J.C. Wolff, J. M. Bioucas-Dias, and M. A. T. Figueiredo, “Near-infrared hyperspectral unmixing based on a minimum volume criterion for fast and accurate chemometric characterization of counterfeit tablets,” Anal. Chem., vol. 82, no. 4, pp. 1462–1469, 2010. [13] D. D. Lee and H. S. Seung, “Learning the parts of objects by non-negative matrix factorization,” Nature, vol. 401, no. 6755, pp. 788–791, Oct. 1999. [14] N. Dobigeon, S. Moussaoui, M. Coulon, J.-Y. Tourneret, and A. Hero, “Joint Bayesian endmember extraction and linear unmixing for hyperspectral imagery,” IEEE Trans. Signal Process., vol. 57, no. 11, pp. 4355–4368, Nov. 2009. [15] O. Gutierrez-Navarro, E. R. Arce-Santana, D. Campos-Delgado, M. O. Mendez, and J. Jo, “A new method to estimate abundances of multiple components using multi-spectral fluorescence lifetime imaging microscopy,” in Proc. IEEE Annu. Int. Conf. Eng. Med. Biol. Soc., 2012, pp. 1081–1084. [16] O. Gutierrez-Navarro, D. Campos-Delgado, E. Arce-Santana, M. Mendez, and J. Jo, “A fully constrained optimization method for time-resolved multispectral fluorescence lifetime imaging microscopy data unmixing,” IEEE Trans. Biomed. Eng., vol. 60, no. 6, pp. 1711–1720, Jun. 2013. [17] C. D. Geddes and J. R. Lakowicz, Topics in Fluorescence Spectroscopy: Glucose Sensing. vol. 11, New York, NY, USA: Springer-Verlag, 2006. [18] A. Periasamy and R. M. Clegg, FLIM Microscopy in Biology and Medicine. London, U.K.: Chapman and Hall, 2009. [19] D. Manolakis and G. Shaw, “Detection algorithms for hyperspectral imaging applications,” IEEE Signal Process. Mag., vol. 19, no. 1, pp. 29–43, Jan. 2002. [20] T.-H. Chan, C.-Y. Chi, Y.-M. Huang, and W.-K. Ma, “A convex analysisbased minimum-volume enclosing simplex algorithm for hyperspectral unmixing,” IEEE Trans. Signal Process., vol. 57, no. 11, pp. 4418–4432, Nov. 2009. [21] P. Pande, B. E. Applegate, and J. A. Jo, “Application of non-negative matrix factorization to multispectral flim data analysis,” Biomed. Opt. Exp., vol. 3, no. 9, pp. 2244–2262, Sep. 2012. [22] F. W. Young, J. de Leeuw, and Y. Takane, “Regression with qualitative and quantitative variables: An alternating least squares method with optimal features,” Psychometrika, vol. 41, no. 4, pp. 505–526, 1976. [23] L. Miao and H. Qi, “Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization,” IEEE Trans. Geosci. Remote Sens., vol. 45, no. 3, pp. 765–777, Mar. 2007. [24] J. Nocedal and S. Wright, Numerical Optimization.. New York, NY, USA: Springer-Verlag, 1999. [25] R. A. Horn and C. Johnson, Matrix Analysis.. Cambridge, U.K.: Cambridge Univ. Press, 1985. [26] J. Jaumot, R. Gargallo, A. de Juan, and Roma Tauler, “A graphical userfriendly interface for MCR-ALS: A new tool for multivariate curve resolution in matlab,” Chemometr. Intell. Lab. Syst., vol. 76, no. 1, pp. 101–110, 2005. [27] B. K. C. Lee, J. Siegel, S. E. D. Webb, L. S. Fort, M. J. Cole, R. Jones, K. Dowling, M. J. Lever, and P. M. W. French, “Application of the stretched exponential function to fluorescence lifetime imaging,” J. Biophys., vol. 81, no. 3, pp. 1265–1274, Sep. 2001. [28] P. Thomas, P. Pande, F. Clubb, J. Adame, and J. A. Jo, “Biochemical imaging of human atherosclerotic plaques with fluorescence lifetime angioscopy,” Photochem. Photobiol., vol. 86, no. 3, pp. 727–731, 2010.

617

Omar Gutierrez-Navarro was born in Tepic, Mexico, in 1984. He received the B.E. degree in electronics engineering from the Universidad Autonoma de San Luis Potosi (UASLP), San Luis Potosi,´ı Mexico, and the M. Sc. degree in computer science and industrial mathematics from the Centro de Investigacion en Matematicas (CIMAT), Guanajuato, Mexico, in 2007. He is currently working toward the Ph.D. Degree in electronics engineering at the UASLP on the study of numerical methods for characterization of living tissue using time-resolved fluorescence lifetime imaging microscopy data. Daniel U. Campos-Delgado received the B.S. degree in electronics engineering from the Autonomous University of San Luis Potosi, San Luis Potosi, Mexico, in 1996, the M.S.E.E. and Ph.D. degrees from Louisiana State University (LSU), Baton Rouge, LA, USA, in 1999 and 2001, respectively. In 2001, he joined the College of Sciences of the Autonomous University of San Luis Potosi as a Professor. He has published more than 120 referred papers in scientific journals and congresses. His research interests include power control in wireless systems, optimization, dynamic modeling, and optimal signal processing. Dr. Campos-Delgado is currently the member of Mexican Academy of Sciences. In May 2001, the College of Engineering of LSU, granted him the “Exemplary Dissertation Award”, and in 2009, he received the “University Award for Technological and Scientific Research” as a Young Researcher from the Autonomous University of San Luis Potosi. Edgar R. Arce-Santana received the B.Sc. degree in computer science engineering from the Technological Institute of San Luis Potosi, San Luis Potosi, Mexico, in 1987) and the M.Sc. and Ph.D. degrees both from the Center for Research in Mathematics in Guanajuato, Guanajuato, Mexico, in 2000 and 2004, respectively. He is currently a Research Professor in the Department of Electronics at the Faculty of Science at the Universidad Autonoma de San Luis Potosi, San Luis Potosi, Mexico. His current research interests include computer vision, optimization, signal processing, and pattern recognition. Martin O. Mendez received the B.Sc. degree in electronics from the Tecn´ologico de Aguascalientes, Aguascalientes, Mexico, in 2001, the M.Sc. degree in bioengineering from the Universidad Aut´onoma Metropolitana, Mexico City, Mexico, in 2003, and the Ph.D. degree from the Department of Bioengineering, Politecnico di Milano, Milano, Italy, in 2007. He is currently in the Department of Electronics, Universidad Aut´onoma de San Luis Potos´ı, San Luis Potos´ı, Mexico, where he is engaged in the analysis and classification of physiological signals to develop clinical decision-support systems. Javier A. Jo received the B.S. degree in electrical engineering from the Pontificia Universidad Catolica del Peru, Lima, Peru, in 1996. He received the M.S. degree in electrical engineering (signal and image processing) and the Ph.D. degree in biomedical engineering (physiological modeling) both from the University of Southern California, LA, USA, in 2000 and 2002, respectively. During 2002–2005, he was a Postdoctoral Fellow in biophotonics within the Department of Surgery, Cedars-Sinai Medical Center, LA, USA. Prior to joining Texas A&M University, he spent a year as a Project Scientist in biomedical engineering at the University of California, Davis, CA, USA. His primary teaching and research interests include systems analysis, signal and image processing, and biomedical instrumentation, with applications to biophotonics, physiology, and medicine.