Introduction

Motion and change are inherent properties of living systems. The time-scales for quantitatively observing dynamic samples can span microseconds for chemical reactions to milliseconds to seconds for cell and organism motion. Technological advances continue to provide access to new regimes in high-speed imaging capable of capturing full volume renderings at rates previously reserved for single frame acquisitions. High speed cameras with kHz frame rates capable of low-light detection have helped propel new methods for biological microscopy [1–4]. High frame-rate 3D sectioning capabilities have also been achieved using temporal focusing for two-photon fluorescence [5] and light-sheet microscopy for conventional fluorescence [6].

Despite these successes, current approaches for high frame-rate imaging using camera-based platforms still suffer from several practical limitations. Firstly, imaging with high frame rates requires high signal to noise ratio (S/N) within each pixel. When detecting sample fluorescence, the high turn-over rates required to produce such S/N can potentially result in significant photobleaching and/or phototoxicity from undesired photochemical reactivity and local heating within the sample [7–10]. In addition, camera-based imaging is not generally compatible with imaging methods that scale nonlinearly with the incident intensity, including two-photon excited fluorescence (TPEF), second harmonic generation (SHG), coherent anti-Stokes Raman microscopy, stimulated Raman gain/loss imaging etc., all of which typically benefit from the higher intensities encountered in beam-scanning instruments. However, one of the most significant limitations of camera-based approaches for high speed biological imaging is the practical difficulties associated with extension to multi-channel detection. Multi-channel detection underpins colocalization experiments, fluorescence resonance energy transfer (FRET) imaging, depolarization ratio detection, and spectral detection. Each channel of detection requires a dedicated high-speed, high sensitivity camera when using the most common high speed microscopy approaches. In addition to the increased complexity and cost of multi-camera detection, precise spatial registry can be challenging to establish and maintain between multiple cameras. In this respect, beam-scanning instruments with multiple single-channel detectors running in parallel offer distinct advantages.

Fast beam-scanning approaches capable of easily supporting multi-channel detection are now well established for video-rate microscopy [11–14]. Beam-scanning is most commonly achieved by combining a slow-scan galvanometer mirror with a resonant mirror or a rotating polygon mirror for video-rate imaging [15]. However, achieving kHz frame rates using a beam-scanning microscope remains challenging. A 512 × 512 image contains ~260,000 pixels, leaving <4 ns per pixel for a 1 ms acquisition to perform the beam positioning and data acquisition. The short duration imposes significant constraints on the beam positioning hardware, the data acquisition electronics, and the sample (e.g., through effects such as saturation, phototoxicity, multiphoton absorption, etc.).

In the present work, a Lissajous trajectory beam scanning technique is implemented on a beam-scanning microscope coupled with a 3D (2D in space, 1D in time) model-based image reconstruction (MBIR) algorithm to recover frame rates exceeding 1000 fps on multiple simultaneous data acquisition channels. By using two resonant fast-scan mirrors with high quality factors (Q > 250) with active phase stabilization for performing the beam-scanning, the position of the beam can be known with high precision at each time-point in the trajectory. Digitizing the signals on multiple channels in synchronicity with the laser repetition rate (80 MHz) allows for reconstruction of high resolution images (> 1 MPix). The MBIR algorithm to recover information in the unsampled pixels builds on the high redundancy in conventional image sets, such that the inherent information content can typically be expressed in much lower dimensionality basis sets.

Adjoining binomial photon counting and photon averaging: theoretical foundation

The underlying Poisson distribution in any pixel can be described entirely by its average value, λ, a good estimate of which is the value that is sought to describe the spectrochemical measurement. The probability of n photons being generated is:

Lissajous trajectory microscopy has been explored previously by several groups in slower imaging applications to achieve frame-rates on the order of a few Hz. Lissajous trajectory beam scanning TPEF microscopes were previously built with MEMS mirrors [16] or fiber scanners [17–19]. Lissajous trajectories have also been utilized in other applications such as atomic force microscopy (AFM) [20, 21]. In general, these previous strategies utilized relatively low least common multipliers or non-repeating Lissajous trajectories with correspondingly low resolution images.

In contrast to other Lissajous studies, the trajectories investigated in the present study utilize several thousand consecutive periods of a resonant mirror prior to the repeating of the trajectory, with a frame-rate that is user definable. Lissajous trajectory microscopy has the distinct advantage over traditional raster scanning systems of allowing simultaneous acquisition of high resolution images at low frame-rates, as well as low resolution images at high frame-rates. The sample is coarsely imaged with only a few passes of the beam, with subsequent passes consecutively adding additional spatial information. Thus, a full Lissajous trajectory image can be temporally divided into several subframes to increase the effective imaging frame-rate with a corresponding tradeoff of in spatial resolution. Figure 1(c) plots the frame-rate vs. pixel fill rate for the 10600:11972 Lissajous trajectory described in this manuscript.

 

Fig. 1 A) Example Lissajous trajectory for 7:8 pattern with equal amplitude. Starting point is 0,0, when the trajectory returns to 0,0 the pattern repeats. B) Pixel sampling density map of 7.083:8 Lissajous trajectory used in imaging. C) Plots of % fill rate vs. frame-rate for the 5300:5986 Lissajous trajectory generating a 512 × 512 pixel image. Higher frame-rates come with the tradeoff of fill rate and/or spatial binning resolution.

Download Full Size | PDF

Model-based iterative image reconstruction

The Lissajous-scanning pattern can leave large blocks of pixels unsampled in each frame at the higher frame rates, although those missing pixels are sampled at other time-points in the trajectory. For kHz regime frame-rates at a 256x256 resolution, less than 20% of the pixels in some frames are sampled. In dynamic systems, the sample can evolve over time, but the signal intensities in the unsampled pixels as a function of position are still correlated with the information in the preceding and following frames. In this way, information acquired as a function of both time and space can be used to more reliably estimate the appropriate intensities for the unsampled pixels, even for objects changing quickly in time. The case at hand can be posed as an in-painting problem as described by Bertalmio et al. [22] Inpainting problems have been solved using a variety of methods, including, but not limited to, partial differential equations (PDEs) [23], patch-based interpolation [24], and sparse learned representations [25], which is an extension to singular value decomposition based on the K-means (K-SVD) [26].

The algorithm used in the current work is unique in that a maximum a posteriori (MAP) estimation with a generalized Gaussian Markov random field (GGMRF) prior model was used, adopted within a 3 dimensional space-time neighborhood. In brief, the MAP estimate identifies the recovered image based on maximizing the probability of observing the measured data given additional imposed model-based constraints. These constraints describe relationships between pixels within a neighborhood. For example, in the simplest implementation of a Gaussian Markov random field, it is assumed that adjacent pixels are normally distributed with a similar mean. However, this simple approach typically does not recover sharp edge effects, resulting in significant blurring in the recovered images. The generalized Gaussian function that serves as the prior model in the GGMRF approach used in this work introduces additional shape parameters in the base Gaussian function, reducing the cost associated with high contrast edges in the images.³² In brief, the exponent in a standard Gaussian function is allowed to vary, with p = 2 (quadratic case) corresponding to a standard Gaussian distribution and p = 1 (total variation case) to an exponential (Laplacian) distribution. It is interesting to note that the GGMRF model in MAP estimation has been used successfully previously on many occasions for recovering sharp edges in low-dose transmission tomography [27] and for blur and noise corrupted images [28], but not in this context of 3D space-time inpainting demonstrated herein.

The developed algorithm is an example of iterative model-based image reconstruction. Intuitively, the prior provides an estimate for the overall probability of obtaining an original image or image set. Next, the algorithm iterates to adjust the estimate for the image, with the degree of change controlled by the regularization. Mathematically, the inverse problem can be succinctly expressed as estimating $\widehat{x}$ from y = Ax + W, where $\widehat{x}$ is the MAP estimate of the original image, y is the observation, A is the Lissajous-scanning matrix, x is the underlying noise free and complete image, and W is noise. While formally the noise in the nonlinear optical measurements is likely to be Poisson-distributed, we will model W as additive white Gaussian noise (AWGN) with a zero mean and a constant standard deviation to simplify subsequent analyses. For high means, the Poisson distribution converges to a Gaussian distribution. Bayes theorem in Eq. (2) provides a starting point for relating the probability of observing the measured results given constraints imposed by the model.

In practice, it is more convenient to minimize the negative logarithm. In this framework, the cost function c(x) to be minimized can be written as shown in Eq. (3),

The first term of c(x) corresponds to the forward model, with the next two terms corresponding to two separate regularizers: one for 2D space and one for 1D time. The relative weighting and scaling of each of the regularizers can be independently adjusted, as they represent distinct physical properties. The scale parameter ${\sigma }_x^p$ can be computed using the maximum likelihood (ML) estimator, which in the asymptotic limit is normally distributed, unbiased, and computationally efficient.

In addition to the iterative model-based image reconstruction approach, analysis was also performed using a complementary interpolation algorithm based on a discrete cosine transform (DCT) coupled with a penalized least-squares (LS) calculation, optimized for implementation in MatLab [29, 30]. A discussion of the DCT-LS algorithm and a comparison between the MBIR and DCT-LS results are included in the Supplementary Information.

Experimental

The Lissajous microscope utilized a Tsunami Ti:Sapphire oscillator (Spectra Physics) pumped by a Millenia Vs (Spectra Physics) CW visible laser operating at 5.0 W. The Tsunami produced femtosecond pulses centered around 800 nm, with an 80 MHz repetition rate, and a maximum power of ~1.0 W. The period of the Tsunami served as the master clock dictating all the subsequent timing of the resonant mirrors, feedback controls, and data acquisition electronics. The incident beam was coupled into the microscope by the first resonant mirror in the resonant mirror pair (Electro-Optical Products Corp.) operating at 5986 clock ticks (13.36 kHz,), providing rapid scanning of the laser beam in the X direction. Using a telecentric lens pair (Thorlabs f = 100 mm each), the incident beam was then coupled to the second resonant mirror operating with a period of 5300 clock ticks (15.09 kHz) to provide fast scanning in the Y direction. The beam was then passed through a second telecentric lens pair providing a 2.3 × beam expansion (Thorlabs f = 25.4 mm, f = 60.0 mm) onto the back of a 10 × objective (Optem, NA = 0.30). The transmitted fundamental beam and the SHG signal were collected and recollimated by a condenser lens (Thorlabs, f = 25.4 mm). The SHG signal was then split from the excitation beam by reflection off a 405 nm long pass dichroic mirror (Chroma, Z405RDC) and passed through a polarizing beam splitting cube (Thorlabs, PBS101). The two polarizations were collected by separate photomultiplier tubes (Hamamatsu, H10722-10) after each passed through filter stacks containing both a 400 nm interference filter (CVI, F40-400.0-4-1.00) and a KG3 (Thorlabs, FGS900) to reject the remaining incident 800 nm light. The transmitted 800 nm fundamental beam was passed through a Glan polarizer and focused by a lens (Thorlabs, f = 25.4 mm) onto a photodiode (Thorlabs, DET10A) for detection of laser transmittance and birefringence. The polarization of the incident light was modulated by an electro-optic modulator (Conoptics, 350-160) placed in the beam path prior to the resonant mirror pair.

Precise timing control was maintained to enable the generation of high resolution images. Due to the high Q-factor of the resonant mirrors (Q > 250) amplitude stability was achieved at the sacrifice of phase stability. To stabilize phase drift, a custom built Lissajous timing generator (LTG) control box was designed to perform real-time active phase correction of the resonant mirrors. The LTG was controlled by an 8-bit Microcontroller (Silicon Laboratories, C8051f120), running at 80 MHz derived from an external 10 MHz phase-lock loop (PLL) synchronous with the 80 MHz master clock from the Ti:Sapphire laser. The microcontroller runs a custom built, multitasking preemptive operating system using a combination of hardware and software timers to produce the X and Y mirror drive signals and epoch pulse, which marks the start of the Lissajous trajectory. The X,Y and epoch signals of the LTG were produced using the synchronous counters of the C8051f120 providing inherent phase lock among these 3 signals. Circuits external to the C8051f120 convert the sinusoidal feedback of the resonant mirrors into a digital signal used by the Microcontroller to compute and compensate for X and Y phase error upon the next iteration of the epoch pulse. User communication with the Microcontroller, to adjust mirror and epoch frequency, was performed via RS-232 with a custom LabVIEW application. The resonant mirrors used in the prototype Lissajous microscope were driven at 5300 and 5986 clock ticks, corresponding to 15.09 kHz and 13.36 kHz, respectively, for an approximate 7:8 trajectory (7.083:8) and a frame rate of ~5 Hz (Lissajous period of ~0.2 s).

Both data acquisition and a significant portion of the timing control was performed using PCI Express digital oscilloscope cards (AlazarTech, ATS 9462), which allowed for continuous streaming of every individual detection event (laser firing) simultaneously on up to 4 channels, although only 3 channels were used to perform the experiments described herein. The digitizer cards were clocked directly from the laser by a fast photodiode (Thorlabs, DET10A), which sampled a small portion of the laser resulting in a pulse train at the frequency of the laser which was used as a clock signal for the digitizer cards. More details about synchronous digitization with beam scanning instrumentation can be found elsewhere [31–34]. The digitizer cards AUX I/O port was configured to supply an output of the clock frequency divided by 8 ( = 10 MHz) in order to clock the LTG through a 10 MHz PLL. As a result the digitization electronics and the mirror drivers all run synchronously with the laser as the master clock (timing diagram, Fig. 2). In this manner, the position within the focal plane, of each pulse of the laser, could be determined with reasonably high confidence. The LTG epoch pulse was used to trigger the start of data acquisition with the digitizing oscilloscope cards.

 

Fig. 2 Timing diagraming for the prototype Lissajous microscope. An 80MHz Ti:Sapphire laser was used as master clock for the ATS 9462 digitizing cards, which in turn generated a 10MHz PLL to synchronize the Lissajous Timing Generator (LTG), which drove the resonant mirror(s) with active phase control. The LTG supplied an epoch pulse to trigger the start of data acquisition.

Download Full Size | PDF

The peak voltages from each detector following each individual laser pulse were digitized and the image(s) were reconstructed according to Eq. (1). Two separate analyses were performed on each single time-dependent data set. In one, each entire trajectory was binned into a single final high resolution image. Alternatively (and simultaneously), the images were separated into sub-frames of the trajectory resulting in higher frame rates of more sparsely sampled images. In either case, the sine-wave trajectory of the resonant mirrors results in non-uniform pixel density across the sample, with higher density near the edges and corners of the image where the resonant mirrors were moving the slowest (turning points). Consequently, some pixels were sampled many times, some fewer times, and in the case of fast sub-frames, some pixels not sampled at all.

Results and discussion

Two model classes of time-varying samples were designed and used for characterization of the prototype Lissajous beam-scanning microscope, selected to represent two distinct effects commonly necessitating high-speed imaging. Sample movement related to fast sample translation (e.g., from sample vibration, heart-beating / breathing, etc.) was modeled by actively translating a sample during imaging. Transient changes in intensity within a fixed field of view (e.g., action potential imaging, vesicle budding, crystal nucleation, etc.) were modeled by performing fast polarization modulation of an SHG-active sample.

In the first model, shown in Fig. 3, the sample was rapidly translated during a 25 Hz Lissajous period image acquisition, such that the spatial position of the sample changed significantly during the course of a single Lissajous trajectory. Figure 3(a) consists of a high resolution image of a USAF 1951 resolution test chart measured by laser transmittance, with no sample translation. Figure 3(b) is of the test chart with rapid translation, resulting in motion blur (Media 1). Figure 3(c) contains an equivalent motion distortion image from traditional raster scanning at 15Hz. The bottom two rows of Fig. 3 explore higher frame-rates with complimentary interpolation as a means of removing motion blur, with 3D and 3F respectively showing sub-trajectories to 100 Hz and 1.25 KHz, and 3E and 3G showing the respective interpolated results (Media 2).

 

Fig. 3 USAF 1951 Resolution Chart. A) High resolution Lissajous microscope image (25 Hz Imaging) of test chart without translation. B) During rapid sample translation, significant blurring occurs. Media 2. C) Example of blurring with traditional raster imaging (15 Hz). D) Data rebinned into sub-trajectories at 100Hz to increase frame rate and reduce blurring with the resultant unsampled pixels. E) Inpainted results of D. F) Data rebinned into sub-trajectories to give 1.25 KHz imaging and the resultant unsampled pixels. G) Inpainted results of F (Media 2).

Download Full Size | PDF

A second model system was developed in which the sample was static in spatial locations, but exhibited transient localized differences in brightness. The measurements, shown in Fig. 4, were achieved by rapid polarization modulation SHG imaging of urea crystals. The period of the polarization modulation was equal to the Lissajous period, resulting in a time-averaged intensity of the field of view for images that are not divided into sub-trajectories. The SHG signal was separated into its horizontal and vertical components with a polarizing beam splitting cube before being detected by fast photomultiplier tubes (PMTs), for sensitive low-light detection (columns in Figs. 4(b) and 4(c)), and complemented by the laser transmittance measurements (column A in Fig. 4(a)) (Media 3–5). The first row shows the full 5 Hz Lissajous period, and the last two rows show sub-trajectories to 1.460 KHz and interpolated results, respectively. The image for the last row was reconstructed based on the sparse sampling provided in the corresponding Lissajous trajectory in the middle row, combined with information from preceding and following spare samplings. The first and last rows are visually similar, though the first row integrates the intensity across the entire modulation period, and the last row features the intensity at a narrow time slice in the modulation period. A movie featuring the 1.460 KHz temporal resolution is provided in Media 3–5, several frames of which are highlighted in Fig. 5 showing the change in intensity as a function of time for the different input polarizations.

 

Fig. 4 Urea Crystals imaged on three different detectors with the incident optical polarization modulated at the Lissajous Period. Columns B and C corresponds to H and V SHG, and column A corresponds to laser transmittance. The first row demonstrates imaging at the 5Hz Lissajous period, the second and third rows demonstrate a sub-trajectory with a 1.460 KHz frame rates and in in-painted results (Media 3, Media 4, Media 5). While the base image in the first row can only show the average intensity across the modulation period, the 292 × increased temporal resolution is able to show the intensity of the various crystal domains at specific time slices.

Download Full Size | PDF

 

Fig. 5 Selected frames of the movie (supplemental information) corresponding image in the last row of Fig. 3 column (A). The enhanced 250 Hz temporal resolution is able to show the intensity evolution of the various crystal domains.

Download Full Size | PDF

The quantitative accuracy of the inpainting image interpolation is most easily assessed in a situation where the ground truth is known a priori. A representative ground-truth video was captured with a CCD camera containing large objects moving in the field of view. Pixels were then selectively deleted from each frame in a simulated Lissajous trajectory for a 250 Hz frame-rate image and interpolated (Fig. 6). For pixels that take on integer grayscale values of between 0 and 255, the interpolated pixels had an average error of 3.4 grayscale units (i.e. 1.5% standard deviation in grayscale value) relative to the ground truth across the interpolated pixels of the image. Comparison with other inpainting algorithms are detailed in the supplementary information.

 

Fig. 6 Movie of “Marco” with toy (A) as a surrogate movie for inpainting interpolation fidelity assessment. After superimposing the same unsampled pixel pattern (B) as in row 2 of Fig. 3, the resulting interpolation is demonstrated in C. For pixels that take on integer grayscale values between 0 and 255, the interpolated pixels had an average error of 3.4 grayscale units ( ± 1.5% standard deviation) relative to the ground truth as displayed in the difference map (D), where gray corresponds to correct values and white and black correspond to over and underestimated vales, respectively.

Download Full Size | PDF

Conclusions

Lissajous trajectory beam-scanning microscopy was coupled with image interpolation algorithms to enable > kHz frame rate imaging on multiple simultaneous data acquisition channels through relatively simple modification of a conventional beam-scanning microscope. Active phase stabilization of the resonant mirrors enabled the acquisition of high-resolution images from Lissajous trajectories several thousand oscillations of the resonant mirrors. By dicing the time-dependent data into smaller sub-trajectory sections, frame rates much higher than those achievable from the overall trajectory could be achieved. MBIR interpolation algorithms enabled reconstruction of the data missing from the unsampled pixels in the sub-trajectory images to produce images at several hundred Hz up to several kHz.

Lissajous trajectory microscopy is particularly advantageous over slow raster scanning techniques in multi-photon and nonlinear optical imaging by minimizing the laser dwell time per pixel. The two most well-established damage mechanisms are due to local heating effects and to undesired photochemistry arising from multi-photon absorption.^{15, 31} Local heating arises from a competition between the rate of heat deposition from the incident beam(s) versus heat dissipation. By using high frequency resonant mirrors in the Lissajous microscope, only 2-3 consecutive laser pulses fall within the same focal volume in the central field of view, with long periods (tens to hundreds of μs) before that same pixel is sampled again in the trajectory. In contrast, conventional imaging using two galvanometer mirrors with 1 μs integration time per pixel and an 80 MHz source corresponds to ~80 consecutive laser pulses, for a >30 times higher local heating rate compared to Lissajous trajectory scanning. In the corners of the Lissajous trajectory, the scan pattern does suffer from an increased consecutive pixel dwell time, where both resonant mirrors are simultaneously moving relatively slowly (Fig. 1(b)). However, adding an iris in the rear conjugate focal plane between the downstream telecentric lens pair provides a simple means to clip just the corners and allow for continuous use of high laser power.

While demonstrated in this study using a pulsed laser source with synchronous digitization, many of the advantages of Lissajous trajectory imaging remain for continuous laser sources (e.g., in confocal fluorescence and confocal reflectance imaging). Given the relative ease with which the experimental hardware can be modified to enable Lissajous microscopy, it should be relatively straightforward to adapt existing confocal microscopes for high frame rate acquisitions by modifications to the scan head and data acquisition electronics.