Analyzing large fluorescence lifetime imaging (FLIM) data is becoming overwhelming; the latest FLIM systems easily produce massive amounts of data, making an efficient analysis more challenging than ever. In this paper we propose the combination of a custom-fit variable projection method, with a Laguerre expansion based deconvolution, to analyze bi-exponential data obtained from time-domain FLIM systems. Unlike nonlinear least squares methods, which require a suitable initial guess from an experienced researcher, the new method is free from manual interventions and hence can support automated analysis. Monte Carlo simulations are carried out on synthesized FLIM data to demonstrate the performance compared to other approaches. The performance is also illustrated on real-life FLIM data obtained from the study of autofluorescence of daisy pollen and the endocytosis of gold nanorods (GNRs) in living cells. In the latter, the fluorescence lifetimes of the GNRs are much shorter than the full width at half maximum of the instrument response function. Overall, our proposed method contains simple steps and shows great promise in realising automated FLIM analysis of large data sets.