Translator Disclaimer
27 February 2019 Adaptive optics microspectrometer for cross-correlation measurement of microfluidic flows
Author Affiliations +
Mapping flows in vivo is essential for the investigation of cardiovascular pathologies in animal models. The limitation of optical-based methods, such as space-time cross correlation, is the scattering of light by the connective and fat components and the direct wave front distortion by large inhomogeneities in the tissue. Nonlinear excitation of the sample fluorescence helps us by reducing light scattering in excitation. However, there is still a limitation on the signal-background due to the wave front distortion. We develop a diffractive optical microscope based on a single spatial light modulator (SLM) with no movable parts. We combine the correction of wave front distortions to the cross-correlation analysis of the flow dynamics. We use the SLM to shine arbitrary patterns of spots on the sample, to correct their optical aberrations, to shift the aberration corrected spot array on the sample for the collection of fluorescence images, and to measure flow velocities from the cross-correlation functions computed between couples of spots. The setup and the algorithms are tested on various microfluidic devices. By applying the adaptive optics correction algorithm, it is possible to increase up to 5 times the signal-to-background ratio and to reduce approximately of the same ratio the uncertainty of the flow speed measurement. By working on grids of spots, we can correct different aberrations in different portions of the field of view, a feature that allows for anisoplanatic aberrations correction. Finally, being more efficient in the excitation, we increase the accuracy of the speed measurement by employing a larger number of spots in the grid despite the fact that the two-photon excitation efficiency scales as the fourth power of this number: we achieve a twofold decrease of the uncertainty and a threefold increase of the accuracy in the evaluation of the flow speed.



The measurement of flows in microfluidic systems or in blood vessels can be obtained either by means of single particle tracking13 or by applying cross-correlation methods on images. Both techniques rely on high-resolution and high-contrast imaging of the sample. Nonlinear excitation methods, in which tightly focused near-infrared pulsed lasers induce the fluorescence emission, have been recently introduced to reduce light scattering from the living tissues4 and increase in general the signal/noise ratio. However, two major phenomena still limit the signal/noise ratio in nonlinear microscopy: the autofluorescence of the tissues and the aberrations of the excitation wave front induced by its propagation through the highly inhomogeneous tissue.

Autofluorescence and scattering that increase the background in different amounts depending on the tissue conservation state4,5 can be reduced by increasing the excitation wavelength. The propagation of the light wave through the tissue increases the pulse width and warps the wave front shape. Both phenomena lead to a decrease of the two-photon excitation efficiency. The pulse width can be reduced by means of a prechirping unit. The wave front distortions can be partially corrected by adding a phase mask on the excitation wave front that precompensates for the tissue induced dephasing within the general approach of adaptive optics (AO).6

Our aim here is to exploit AOs methods for the correction of the optical aberrations in order to improve cross-correlation analysis of flows in biological samples. Cross-correlation methods, in which we measure the time-of-flight of fluorescent microbeads between observation volumes lying along the flow lines, require parallel (i.e., simultaneous) collection of the signals from sparse observation volumes with high pointing stability and high signal/background ratio. Our rationale is to use a single diffractive programmable element to obtain multispot excitation of the sample and to correct optical aberrations on the spots in order to improve the quality of the interspots cross-correlation function (CCF). Our results depend on the unique coupling of cross-correlation image analysis and aberration correction, achieved on a single optical path without any movable mechanical part. Therefore, we briefly review these two major research areas in the following.

Spatial light modulators (SLMs) have been widely used to change the image modality7 of optical microscopes as well as to correct aberrations. For example, Gharbi et al.7 implemented a simple method based on multiple carrier frequencies to reduce the chromatic aberrations in an incoherent illumination microscope. Hasler et al.8 developed a two-image recording method with digital postprocessing based on an SLM, which is equivalent to a phase contrast microscope. Lingel et al.9 developed an optical setup for phase retrieval, in which the object can be masked first in the intermediate image plane and then filtered in the Fourier plane by making use of a single SLM and no movable mechanical parts. Even more recently, Mastusmoto et al.10 described a simple algorithm to correct the large aberrations induced by the shape of the surface of the biological tissues in two-photon excitation microscopy with a dry objective.

The possibility to implement AOs corrections in raster scanning setups for nonlinear microscopy has been first reported by Neil et al.11 These authors presented an intensity-based sensor-less AOs method for two-photon raster scanning microscopy obtained by implementing a differential mapping of the phase of the object. Débarre et al.12 substantially improved on this scheme and demonstrated the feasibility of SLM-based aberration correction of two-photon images of biological specimens.

The deformable mirror and multiactuator deformable lens are possible alternatives to the SLM. Two general approaches to the wave front correction are reported in the literature: zonal correction and modal correction13,14. In the first case, the wave front is described as the superposition of Ξ basis functions, typically the Zernike polynomials, and the correction is achieved by changing their amplitudes. In the second case, the correction is performed on the back pupil plane of the objective, which is divided in S regions of interest. Deformable mirrors are not much sensitive to the light polarization and wavelength15 and allow for modal or zonal correction of the aberrations.13,14 Multiactuator deformable lenses easily allow for modal correction and depend on the light wavelength.16 The SLM offers a huge number of degrees of freedom for the correction of the aberrations at the cost of a slow refresh rate. It allows to apply modal or zonal corrections of the aberrations and to create customized illumination patterns, such as grids of spots, on which cross-correlation analysis can be easily performed.

Cross-correlation analysis of flows in microfluidic devices constitutes a wide research field with implementations that span from microPIV data analysis17 to flow mapping on single images.18 It has recently acquired great impact in cellular biology, where it is applied on 2-D images19,20 in order to characterize the flow field in the image plane.21 The cross-correlation analysis of the sample dynamics can be based on a multispots simultaneous excitation of the sample2225 or on a spatially distributed and time multiplexed excitation imaging as in raster image correlation spectroscopy.8,26 Multifocal excitation and collection, occurring in parallel on a pixelated device, offer best results for the analysis of fast flows. Multifocal excitation was also used to increase the image acquisition rate or to gain a parallel acquisition on a sparse array of observation volumes. The laser beam was mechanically moved on the sample and coupled to an acousto-optic deflector2729 or an SLM30 in order to obtain a multispot excitation of the sample. In these works, a fine adjustment of the position of the spots on the spines was done by means of galvanometric mirrors. However, the pointing stability was not assessed directly or by means of a cross-correlation analyses.

Based on these works, we thought that only scan-less setups with no movable parts could offer the parallel acquisition needed to perform cross-correlation analysis. Some of us preliminary tested this assumption in collaboration with D’Angelo et al.31 on arrays of neurons in acute cerebellar slices of rats. The scan-less parallel acquisition developed in Ref. 31, based on the use of an SLM, offered the possibility to draw connectivity maps of neurological activity in an ensemble of communicating neurons32 in the granular layer of acute rat cerebellar slices. The images were collected by shifting an array of laser spots on the sample, by collecting the signal of the whole array of spots through a sensitive EM-CCD camera, and by superimposing these signals on a single map. The number of spots in the array and the shifting step size determined the frame rate and the image resolution: the larger was the number of spots per array, the lower was the image acquisition time. However, the number of spots per array cannot be raised without rapidly degrading two-photon excitation efficiency, a fact that could be partially overcome by optimizing the laser beam shape. In Ref. 31, the SLM was not exploited to correct aberrations of the optics that could result in an increase of the two-photon excitation efficiency. In addition, no attempt was done to measure the CCFs between the spots of the array.

The novelty of the present report with respect to Ref. 31 is, first, to substantially improve the localization of the spots on the sample (or camera) plane and implement a more accurate biquadratic mapping method between the SLM and the camera planes. This is essential to obtain accurate mapping of flows and in general to reduce artifacts in cross-correlation analysis on sparse samples. Second, we implement on the same SLM an optical aberration correction algorithm that increases the performances of the microspectrometer for combined morphological and fluidodynamics studies. The SLM is used to shine rectangular arrays of spots, to correct their optical aberrations, and to shift them on the sample so to collect two-photon excitation images of biological specimens. In a second step, the SLM is used to shine arbitrary patterns of aberration corrected laser spots, on which we compute the CCFs that provide the measurement of the flow velocity. This setup is, therefore, capable to collect aberrations-corrected two-photon images, to point at specific regions of interest with high accuracy, and to increase the sensitivity of the cross-correlation analysis of flow velocities. We combine for the first time the correction of wave front distortions to an excellent pointing accuracy on arbitrary geometries of excitation volumes, achieving a fivefold increase in the signal-to-noise ratio of the acquired CCFs, a twofold decrease of the uncertainty and a threefold increase of the accuracy in the evaluation of the flow speed.


Materials and Methods


Optical Setup

The optical setup consists of a custom-made upright two-photon holographic microscope (Fig. 1). Two-photon excitation is obtained in the samples by a Ti:Saph laser (Tsunami, Spectra Physics, California) providing 1 W maximum power at λ=800  nm. The polarization of the laser light is perpendicular to the optical table. The laser is collimated and expanded by a 3.5× beam expander (f1=40  mm and f2=140  mm) and impinges at 9 deg on the reflective SLM (X10468-7, Hamamatsu, Japan), which allows to encode phases in the range 0 to 2π on the beam wave front with 792×600  pixel resolution (20-μmpixel size) and 8-bit pixel depth. The SLM plane is conjugated through a 0.65× beam reducer (f3=200  mm and f4=125  mm) to the entrance pupil of a microscope objective (20× XLUMPLFLN, Olympus, Japan). An arbitrary intensity distribution in the focal plane is obtained by encoding a phase mask on the SLM.33,34 The contribution of the zero-order diffraction from the SLM onto the sample plane is minimized by superimposing a Fresnel lens (FL) pattern on the calculated phase masks,35 while slightly displacing the position of the second lens of the beam reducer from the telescope condition.31

Fig. 1

(a) Experimental setup. The lenses L1 and L2 expand the laser beam to fill the SLM. Lenses L3 and L4 act as a beam reducer to fill the entrance pupil of the objective (vertical); L3 and L4 are placed at a distance lower than the sum of their focal length, to allow the elimination of zero-order contribution. The fluorescence is collected through the DM, filtered by a bandpass filter (F), and focused on the CCD by the tube lens (LT). Two folding mirrors (FM) are used in the setup. (b) Schematics of the image acquisition algorithm in the CMOS plane. The symbols refer to Eqs. (9)–(12) in the text. The red square (upper left) in the reconstructed image corresponds to I=J=3 and k=l=0 and i=1, j=1 in Eq. (12).


The sample fluorescence is discriminated from the excitation light through a dichroic mirror (DM) (HC670SP, AnalysenTechnik, D) and a bandpass filter (ET680SP-2P8, AnalysenTechnik, D). An image of the objective focal plane is formed by a LT (fLT=150  mm) on a Hamamatsu CMOS Orca Flash 4.0 Camera (pixel size 6.5  μm).


SLM-Based Illumination Pattern Creation

Aberrations that detrimentally affect the resolution and the contrast of microscope images can be corrected by means of AO with a variety of active diffractive or reflective elements.11,36,37 Here we exploit an SLM to create custom grids to implement the cross-correlation methods on a multiple spots basis and to acquire images with no movable parts in the optical path similar to Ref. 31. Differently from that work, we combine here these features with the possibility offered by the SLM to correct the wave front distortions on different regions of interest throughout the field of view.


Grid Projection

A grid of N×N spots is created in the sample plane by loading onto the SLM a phase pattern computed through the weighted Gerchberg-Saxton (GSW) algorithm.33,34 The resulting phase pattern (GSW-phase) is modulated between 0 to 2π before applying it onto the SLM. The grid intensity profile projected onto the sample plane is rescaled on the two axes, since a discrete 2-D Fourier transform is applied on a rectangular domain (the SLM plane) with different sizes (792×600  pixel). We calculate two conversion factors κcol and κrow, respectively, equal to the ratios between the mean horizontal and vertical spacing between two spots in the SLM domain and in the camera plane. These factors allow calculating the horizontal and vertical spacing in the SLM domain to produce a squared grid in the sample plane [see Eq. (10)].


Zero-Order Suppression

One of the major issues related to the use of the SLM is the so-called zero-order background.35 A percentage (<25%) of the collimated light incident onto the SLM is not modulated by the pixelated structure and produces non-negligible levels of diffused light in the sample plane, which results in a limited signal-to-background ratio of the SLM-generated illumination pattern. In order to reduce the zero-order contribution, we placed the L3 and L4 lenses [Fig. 1(a)] at a distance shorter than the sum of their focal lengths of a quantity |2ε|, ε<0, and we add a quadratic phase pattern (FL) to the GSW-phase.31 The FL has the following expression:

Eq. (1)

where κ=n/λ and f is the focal length of the lens L3 in Fig. 1(a) (n = refraction index = 1). An FL phase-modulated wavefront passing through the lens L3 is focused on a plane placed at a distance 2ε with respect to the lens focal point, but corresponding to the focal point of the lens L4. Therefore, it will leave the beam reducer as a collimated beam. The zero-order contribution will instead diverge since it is focused after the focal point of the lens L4 along the optical path.


Setup Calibration

The image reconstruction algorithm is based on a grid of spots that is generated by the SLM and shifted over the field of view. The cross-correlation method is applied between couples of spots on the sample plane. In both cases, it is fundamental to identify with high accuracy of the position of the spots in the camera plane. Compared to Ref. 31, we improved substantially the accuracy of the localization of the laser spots on the sample by means of blob-detection algorithms (BDAs)38 and of a biquadratic mapping between the SLM and the camera planes. We checked that after a careful calibration of the sample-to-CMOS planes mapping the spots positions on the camera could be identified for a wide variety of grid parameters.


Mapping calibration

We start the calibration algorithm by projecting a grid with known features onto a Rhodamine plastic slab and acquiring the fluorescence signal with the CMOS camera. In order to identify the groups of pixels that correspond to each spot of the grid, we exploited a BDA:38 a blob is a group of connected pixels that share some common properties, e.g., grayscale value. The BDA works as a spatially resolved bandpass filter on the signal levels between a minimum and a maximum threshold: in our case, the maximum threshold was set to 255 digital levels (8-bits depth). It must be noticed that the position of the blob in the BDA is a real variable not limited to the integer values that identify the real pixels. Since the grid was typically affected by an intensity gradient, we first subtract from the grid image the zero-order contribution obtained by deleting the FL mask from the GSW algorithm. The image was then split in several subregions, the BDA-detection algorithm was applied separately on each of them, and the detected spots were then stored in an array. Finally, for each position stored in the array, we copied a regions of interest (ROI) centered at the putative spots from the image of the zero-order subtracted grid image. We then applied a Gaussian filter and autocontrast to the whole new image before iteratively applying the BDA for increasing values of the minimum threshold from 0 to 255 until we found a number of blobs equal to the number of spots of the grid. As a final step, we took as putative spot position in the camera plane the maximum pixel coordinates value inside the ROI centered on the blob position.


Prediction of spots position

In order to estimate the position of the spots of a grid with features different from the calibration ones, we adopted a 2-D second-order parametric model39 that describes the mapping of the SLM plane onto the sample plane. This model allows us to retrieve the spots position without applying the BDA each time we change the grid shape. The transformation between the position of a single spot in the SLM domain (xs,ys) and the camera domain (xc,yc) can be expressed as

Eq. (2)


We adopt a matrix notation, where P is a N×2 matrix (N is the number of spots in the grid) whose rows contains the coordinates (xc,i,yc,i) of the i’th spot in the camera domain, C is a N×9 matrix whose i’th row contains the products {xs,inys,im}m,n=0,1,2, and A is a 9×2 matrix that contains, in each row, the expansion coefficients of xc,i and yc,i into polynomials [Eq. (2)]. Once the matrices P and C are obtained from the mapping calibration procedure, the transformation matrix A is simply calculated as39

Eq. (3)


The matrix A is computed once from the data obtained in the mapping calibration procedure for a specific choice of the grid parameters. For different grid spacings in the sample, characterized by a A matrix, the position of the new spots in the camera domain could be estimated as

Eq. (4)



Aberrations Correction

Our aberrations-correction algorithm is based on a sensorless AO scheme.11,17,40,41 We assume that a wavefront can be expressed as a weighted sequence of basis functions or aberration modes and that the image quality can be restored by maximizing a metric that depends on the expansion coefficients of the wave front in terms of the modes.


Aberrations representational model

The wave front propagating toward the sample plane can be described at the entrance pupil as a weighted sum of basis functions:

Eq. (5)

where ψi is the i’th basis function and zi is the corresponding amplitude. The most common choice over a circular aperture is the Zernike set of polynomials, which is a complete and orthogonal set of functions defined over a unit circle.42 We needed here to extend this basis to the Cartesian representation of the Zernike modes43 due to the rectangular symmetry of the SLM. In order to do so we followed the algorithm described in Ref. 43 and adopted a Noll’s index scheme.44 We have limited the aberration corrections to the 11’th Zernike polynomial (primary spherical).

The versatility of our system in creating custom illumination pattern is challenged by its low-speed performance in acquiring images. For this reason, we decided to build a wave front metric based on the intensity measured at selected spots on the sample plane, instead than on the whole image as done in other studies11,45,46 based on raster scanning or wide field image acquisition. As we discuss in Sec. 3, this algorithm could also be modified to account for the correction of nonisoplanatic aberrations.47,48 Our aim is to obtain bright and small spots on a low zero-order background after the application of the aberrations-correction algorithm.

Our metric is based on the comparison of the intensities Sn^ integrated over ROIs (size m×m) centered on each of the N spots of the grid ({xc,i,  yc,i}i=1,,N) with the background measured by calculating the mean integrated intensity B^ over four ROIs (size m×m) set at the vertexes of a square centered on the grid central spot whose side is equal to the spot–spot distance. Throughout this report m=3, unless otherwise specified. By limiting the computation of the background to the central spot, we are overestimating the effect of the zero-order reflection from the SLM that is the primary source of background and fades at the edges of the field of view. Therefore, we computed the metric over the whole field of view covered by the N spots as

Eq. (6)



Metric optimization

A common assumption, well posed in many practical situation, is to describe the wavefront as a finite sum of Ξ Zernike modes with amplitudes {zi}i=1,,Ξ and to approximate the metric with a paraboloid surface in the vicinity of its maximum:11,40,45,49

Eq. (7)


In Eq. (7), q and γ={γkl}k,l=1S are the constants. This function has a clearly defined maximum, but each term of the sum depends upon the amplitude coefficients of two aberration modes. Since the matrix γ is not strictly diagonal, the optimization of the metric is not straightforward. It is possible to diagonalize the matrix17,49 γ and get, at high computational cost, a correction improvement that depends on the amount of intermode coupling. As proposed in previous works,11,17,45 we adopted two searching algorithms based on the maximum parabolic interpolation.50 The first one searches sequentially the maximum of the metric as a function of the individual k’th mode amplitude, while depressing the effect of the other (zl=0 if kl). This leads to the best mode amplitude (BMA) for each of the considered Zernike polynomial. We then sum the corrections of all the modes at their BMA. In an alternative algorithm, we searched the global maximum of the metric by means of the gradient descent method (GDM). In this algorithm, we initialize all amplitude modes to zero, numerically calculate the gradient of the metric M by means of the Newton’s difference quotient method, exploit the backtracking line search algorithm51 to define the optimum step to update all the amplitude modes, and iterate this procedure until the absolute value of the gradient is lower than an arbitrary chosen cut-off.

Experimentally, the two methods produce similar results in the samples investigated here. This means that the mode that at its BMA produces the highest correction to the metric is typically the mode whose relative derivative at the starting point of GDM is larger. In general, by independently maximizing the metric on each mode, we obtain an estimate of its contribution to the wave front aberrations. Instead, the GDM provides an aberration correction pattern that involves the contribution of all the modes and may lead to a slightly larger improvement of the metric.


Image Reconstruction

Our aim is here only to characterize the morphology of the sample before the flow measurement without any moving optical parts in order to be able to position with high accuracy of the observation volumes for the cross-correlation measurements and to measure the flow velocity. The image is acquired through a shift-and-snap procedure.31 We project a Nrow×Ncol spots grid onto the sample plane. Then the signal of each spot is read and stored in the corresponding pixel of the image [the position of each spot in the camera plane can be predicted by Eq. (4)]. We shift the grid in the sample plane along horizontal and vertical directions and iterate the procedure, building up the final reconstructed image. The grid is shifted along the horizontal direction by adding a prism phase function to the grid pattern profile:

Eq. (8)

that moves the spot of a quantity |δP|λf3 in the focal plane of the lens L3 [Fig. 1(b)]. For the shift along the row direction δP=(δrow,0), while for the shift along the column direction δP=(0,δcol).

The number of shift movements along horizontal and vertical directions (respectively, ΛH and ΛV) is determined by the number of spots and the size in pixels of the reconstructed image:

Eq. (9)

where W and H are the reconstructed image dimensions [width and height, see Fig. 1(b)]. The spacing between spots in the SLM domain ΔcolSLM, ΔrowSLM depends on the number of shifts and the pixel size in the image plane (δcol and δrow):

Eq. (10)

where M is the magnification factor between the sample and the camera planes, κcol and κrow are the grid-distortion correction factors (see Sec. 2.3) and Px, Py are the linear sizes of the camera pixel in μm (in our case Px=Py=6.5  μm). The shift of the grid (in terms of SLM pixels in the SLM plane) can be determined as

Eq. (11)


Finally, the signal read from the spot in the i’th column (from left to right, i=0,1,,Ncol) and j’th row (from top to bottom, j=0,1,,Nrow) of a grid that has been shifted l times along the horizontal direction (l=0,1,,ΛH) and k times along the vertical direction (k=0,1,,ΛV) is assigned in the synthetic reconstructed image to the pixel with indices {I,J}:

Eq. (12)


When we reconstruct an image, we usually integrate the signal from a 3×3 ROI centered in every putative spot position, in order to minimize the fluctuation of the spot intensity. The description of the grid shift algorithm given above in the SLM plane by means of the spot–spot distance (ΔcolSLM,ΔrowSLM) and the shift (δcolSLM,δrowSLM) parameters can be converted in the sample plane by means of equivalent parameters called (Δx,Δy) and (δx,δy) used in Sec. 3.


Cross-Correlation Measurements

The flow velocity can be measured from the analysis of multiple CCFs. The two-volumes normalized CCF is defined as

Eq. (13)

where the signals I1(t) and I2(t) are collected from two observation volumes spaced by a vector R. For colloidal microparticles (or cells) drifting with velocity V, we assume that the speed |V| is constant between the two volumes and that it makes an angle α with the vector R. The CCF is a quasi-Gaussian decay modulated by a hyperbolic decay function:4,6,52

Eq. (14)

GX12(τ;τD,τv)=GDiff  (τ;τD)·exp{|v|2ω02[τ2+τv22ττvcos(α)]1+τ/τD},GDiff(τ;τD)=1N11+τ/τD1+ξτ/τD.

In Eq. (14), τD=ω02/4D is the diffusion time, D is the diffusion coefficient of the colloids, ω0 is the beam waist, ξ1020, accounts for the volume elongation along the optical axis, and τv=|R|/|V| is the time of flight between the observation volumes. When the flow axis is collinear (α=0) to the interspot distance R, Eq. (14) can be approximated by a simple Gaussian function in the lag time with maximum at the lag time τmax=τv, for small values of the times ratio Γ=τv/τD. Numerical simulations proved6 that for Γ<0.1 the discrepancy between τv and τmax is lower than 0.5%. If the interspot vector is not collinear to the flow velocity, the argument of the exponential function in Eq. (14) does vanish exactly, the lagtime of the maximum of the CCF is displaced with respect to τv as τmaxτvcos(α). Moreover, the amplitude GX12(τmax;τD,τv) of the CCF decreases as a function of the noncollinearity angle α as GX12(τmax;τD,τv)exp(α2/β2), where β is the angular size of the downstream volume as seen by the upstream volume. The correction to the position and the amplitude of the maximum of the CCF is second order in the misalignment angle α.

Equation (14) provides only the flow speed, through the measure of τv, in which the maximum of the CCF falls. In order to map the flow velocity, we projected grids of various symmetries on the sample plane and calculated the CCF between each pair of neighboring spots. For each spot (called hereafter C-spot), we retrieved the neighboring spot (called hereafter NN-spot), for which the CCF shows the highest cross-correlation amplitude and assign to the C-spot a velocity whose amplitude is the speed calculated from τv and whose direction lies along the vector pointing from the C-spot to the NN-spot.

In an effort to enhance the signal/noise of the cross-CCFs, we applied clipped correlation algorithm,53 in which the time trace of the signal collected on an ROI is digitalized according to a threshold value, Th. All the time points for which the signal exceeds Th are set to 1, otherwise they are set to 0.




Spots Identification and Spot Prediction Algorithm

Both image reconstruction and cross-correlation analysis of flows require an accurate identification of the spots in the grid. We validated the stability of our spot identification algorithm on a variety of grids with different parameters (spots number and spot–spot spacing) and for different values of the signal-to-noise ratio. An example of the algorithm output is given in Fig. 2(a), for a 5×5 grid (spot–spot spacing 21.1×15.8  μm in the sample plane). In this example, the signal-to-background ratio was 9.5±0.4 and the pixel with the highest intensity coincides within one pixel with the putative spot position obtained by means of the BDA algorithm (see Sec. 2.5.1). By analyzing a series of low signal/noise ratio images obtained by reducing the laser power or the camera exposure time, we checked that the BDA algorithm retrieves the spots position well within one pixel (average positioning error  50  nm17  pixel) for signal/background ratios as low as 1.35 [see Fig. 2(b)]. The accuracy of the spot prediction algorithm was validated also on shifted grids by applying the spot retrieval algorithm on grids projected onto a Rhodamine plastic slab, finding similar accuracy over the whole range of shifting values employed here.

Fig. 2

(a) Results of the spot identification algorithm on a two-photon image of a 5×5 grid projected onto a Rhodamine slab (interspot spacing 21.1×15.8  μm, sample plane). The grid was corrected for optical aberrations. Each identified spot is represented by a red circle centered in the corresponding putative spot position. (b) Root-mean-square error in the assignment (BDA algorithm) of the spot center as a function of the image signal/noise. (c) Demonstration of aberration correction efficacy on the two-photon fluorescence image of a 5×5 grid (interspot spacing 28.8×21.8  μm, sample plane) projected onto a Rhodamine slab (upper-left, uncorrected image; lower-right, aberration corrected image). Both the signal-to-background ratio and the width of each spot improved after aberration correction. The plot profile of the dashed line in (c) for the uncorrected (solid line) and aberration corrected (dashed lines) grids are reported in (d).



Aberrations Correction


Algorithm for correction on a single array of spots

We analyzed and corrected the optical aberrations affecting a grid of spots projected onto a Rhodamine plastic slab by choosing at first to maximize the metric [Eq. (6)] independently on each Zernike mode and tested modes up to n=11 (primary spherical). The first estimate of the aberration corrections is the ratio between the metric of the uncorrected and corrected image. In order to enhance the image quality, we superimposed the corrections for the modes that exhibit the metric improvement 20%. The improvement in the two-photon image of a 5×5 grid (interspots spacing 28.8×21.8  μm) projected onto a Rhodamine slab can be judged visually from Fig. 2(c). The analysis of the images [line profiles reported in Fig. 2(d)] indicates that the relative improvement of the signal-to-background ratio between the raw (signal/background=7.3±0.6) and the corrected spots (signal/background=2.4±0.1) lying on the dashed line in Fig. 2(c) is 3.0±0.3.



We studied the dependence of the aberrations on the position of the grid within the FOV and the possibility to apply different corrections over a wide FOV. This was done by shifting a 5×5 grid in different regions of the FOV onto the Rhodamine slab [Fig. 3(a)]. The correction amplitude of each mode [zn in Eq. (5), 5n11] changes over the ROIs as it is reported in Fig. 3(b). The analysis reported in Fig. 3(a) and 3(b) could be performed for an extended field of view in a thick biological sample, in which scattering from different regions of the out-of-focus portions of the tissue produces nonuniform aberrations, and therefore, requires spatially dependent corrections.

Fig. 3

(a) Juxtaposition of nine 5×5 grids (aberrations corrected) projected onto a Rhodamine slab (two-photon fluorescence, interspot spacing 21.1×15.8  μm in the sample plane) with different combinations of x shift (131.0, 0, and +131.0  μm in the sample plane) and y shift (99.8,0, and +99.8  μm in the sample plane). (b) For each ROI from 1 to 9, the mode amplitude (left y axis of each of the 9 plots) and the corresponding maximum metric value (right y axis) are shown for each Zernike modes from 5th to 11th (x axis). (c)–(e) Two-photon images of 5×5 grids projected onto a Rhodamine slab with increasing interspot spacing: (c) 7.7×6.1  μm, (d) 34.8×26.3  μm, and (e) 67.5×52.7  μm. The inset in (c) reports the profile along the central horizontal line of the grids in (c) (light gray), (d) (dark gray), and (e) (white).



Signal/background and spots size

The signal/background of the grid (average of the maximum of the spot divided by the background between two adjacent spots) is increased typically by a factor of 5±1.5 times by the application of the Zernike modes corrections. However, the efficiency of the aberrations correction is largely affected by the spot–spot spacing in the grid. The signal/background ratio decreases rapidly with the reduction of the spacing among the spots of the grid: by changing the spot–spot spacing from 67 to 8  μm, the signal/noise ratio on the corrected image decreases from 80±10 to 3±0.7 [Figs. 3(c)3(e)].

The effect of the aberration correction on the single spot is a considerable reduction of its full width passing from 4.5±0.8  μm to 2.67±0.02  μm for the same number of spots in the matrix [5×5, Fig. 4(a) middle plot and Fig. 4(b)]. However, the aberrations-correction effect on the width of the spots depends on the number of spots of the grid at fixed spot–spot distance [Fig. 4(a)]. We report as an example in Fig. 4(a) the trend of the spot width as a function of the number of spots for images taken at spot–spot distance=27.8  μm. By changing the grid from a 3×3 to a 7×7 matrix, the full-width at half-maximum decreases from 2.90±0.03 to 2.33±0.02  μm [Fig. 4(a)].

Fig. 4

(a) Single-axis profiles of spots obtained with aberration corrected grids of increasing number of spots 3×3, 5×5, and 7×7, respectively, from bottom to top, each with interspots spacing equal to 27.8×20.7  μm. The full-width at half-maximum of the spots measured in these three cases is 2.90±0.03, 2.75±0.03, and 2.33±0.03  μm, respectively. (b) The radial distribution for a noncorrected 5×5 grid (same interspot spacing) has a full-width at half-maximum=4.5±0.8  μm, almost double the nonaberrated case. (c) and (d) Two-photon images of a Convallaria sample obtained by employing a 13×13 grid (spot–spot spacing 27.4×27.4  μm) before aberrations correction [BC, (c)] and after the aberrations correction [AC, (d)] field of view=356×356  μm, excitation power P=150  mW, acquisition time per frame = 10 ms, total number of frames acquired to reconstruct the image = 1521. (e) Direct comparison of two-photon images of a Convallaria sample before (BC) and after (AC) the aberrations correction. A 5×5 grid has been employed with a spacing Δx=14.0  μm and Δy=14.0  μm in the sample plane, shift steps dx=702  nm and dy=702  nm in the sample plane, corresponding to 2 CCD pixels in the CCD plane for a FOVs of 70.3×70.3  μm: the excitation power was P=150  mW, the acquisition time per frame was 10 ms, and the total number of frames acquired to reconstruct the image was 400. (f) Profiles of the Convallaria image (e) along the dashed and the solid white lines marked in (e).



Image Reconstruction and Correction

We tested our image reconstruction and aberration correction algorithm on a sample of Convallaria (Convallaria majalis stained with fast green safranin, AS3211, Leica, D). By projecting a grid of spots and shifting it over the sample in the x and y directions, we were able to reconstruct the nonlinear excitation autofluorescence image (excitation at λ=800  nm) of the sample. Figure 4(c) shows an example of the images obtained before the aberrations correction, by means of a 13×13 spots grid (spacing Δx=Δy=27.4  μm in the sample plane) shifted with steps δx=δy=702  nm in the sample plane. The grid shift corresponds to 2 CCD pixels in the camera plane. The grid rectangular distortion was quantified and corrected with correction factors κcol=3.0±0.05 and κrow=4.0±0.05 in the x and y directions, respectively. For each pixel of the reconstructed image, the intensity was computed as the mean of the signal on a 3×3 mask centered in the corresponding spot of the grid employed in the specific image reconstruction step. The efficacy of the spot array aberration correction can be visually estimated over the whole field of view from the comparison of Figs. 4(c) to 4(d) and better evaluated from the zoom reported in Fig. 4(e), where uncorrected and aberration-corrected images are juxtaposed on the same image. To correct the optical aberrations affecting our setup, we evaluated the metrics of Zernike modes from 5th to 11th on a Rhodamine slab, optimizing all the modes simultaneously (gradient search GDM, see Sec. 2.6.2) on a fixed 13×13 grid of spots. We then collected the image by summing the aberrations-correction phase pattern to the grid phase pattern at each acquisition step. The mean signal-to-noise ratio increases more than twofold (2.2±0.1), as can be judged from the profiles reported in Fig. 4(f). The resolution of the structures of the sample is clearly better after the aberrations correction as can be judged from the possibility to resolve the y-shaped structure in Fig. 4(e) along the solid and dashed white lines for the aberration corrected (AC) and uncorrected (BC) images at about x=29 and 33  μm, respectively, [Fig. 4(f)]. The full-width at half-maximum of the two peaks that are clearly discernible in the aberration corrected image profile is 2.1±0.1  μm. High-resolution imaging of the same sample features with a Leica SP5 confocal microscope (resolution 330±30  nm) gives a width for these structures of 1.4±0.5  μm. These values indicate an optical resolution of the present setup of 1.5±0.5  μm.


Cross-Correlation on Grid of Spots for Flow Mapping

After getting an image of the sample, the grid array can be positioned at a specific position with an accuracy that is better than a single pixel [see Fig. 2(b)] and any bright object moving through the spots can be detected with a resolution of 10 ms, determined by the camera acquisition time. We tested the accuracy with which we can measure the flow speed on microchannels of varying geometries, in which we injected microbeads suspensions at fixed speed values. The microbeads passing through the spots are detected as bursts of signal on the camera image. The density of the microbeads is low enough to provide distinct sparse bursts of signal on the time profile of the time stack of images.

Two geometries of microchannels were studied: a straight square cross-sectional glass microchannel (800  μm×800  μm size) and a horseshoe-shaped microchannel (500  μm×500  μm size) obtained by PDMS molding. We exploited the full capability of the SLM in creating custom illumination profiles and measured the cross correlation of the signals collected from the two volumes (the signal was integrated on a 3×3  pixels area). The principle of cross-correlation spectroscopy52,54 requires the synchronous collection of the signals from two laser spots along the flow, a possibility offered by the SLM that also allows to maximize the signal/noise ratio by correcting the optical aberrations on the selected spots.


Straight microchannel

At first, we applied the cross-correlation method on couples of spots lying along the flow of a straight glass microchannel. The two-photon excitation depends on the square of the excitation intensity. By multiplexing the laser power on a N×N grid, we reduce the excitation intensity on each spot of the grid by a factor scaling as N4. So, as a first choice, a 5×2 rectangular grid of spots (Δx=14.8  μm and Δy=4.6  μm, sample plane) was projected on the sample. We estimated the direction of the flow by cross correlating the central spot [red dashed circle C’, Fig. 5(a)] with each of its nearest neighbors [orange circles, from A to E, Fig. 5(a)].

Fig. 5

CCF analysis on a portion of a straight microchannel. (a) Image (maximum signal projection of a time stack) of a 5×2 grid (interspot spacing 14.8×4.6  μm, sample plane). (b) Cross correlations between the C’ spot and the A–E spots [as in (a)], respectively, without (top) and with (bottom) the correction of optical aberrations. The solid black lines are the best fit of the data to a Gaussian function for the C’–C couples of spots with the best fit variances of 2.5±0.2 and 1.6±0.1  ms for the noncorrected and aberrations corrected cases, respectively. (c) Increase of the cross-correlation amplitude at the time-of-flight (τmax) as a function of the cut-off threshold (ROI size for the signal average, 3×3  pixels). The inset reports the CCFs for increasing values of the threshold. (d) Effect of the size (m×m) of the ROI used to average the signal on the spot (threshold for cut-off = 25): 1×1 (squares), 3×3 (circles), 5×5 (up-triangles), and 7×7 (down-triangles). The solid lines are the best fit Gaussian functions. The left and right panels report the normalized CCFs as a function of the lag time and the amplitude GX12(τmax) as a function of the ROI size in pixels, respectively. (e) CCFs calculated for the C–C’ pair of spots (a) for different flow speeds (set values=120, 400, and 1200  μm/s), for the case of no correction of the optical aberrations (solid lines), and for the case of aberrations correction (dashed lines). (f) and (g) Maximum projection images of a 4×4 grid (interspots spacing 7.0×5.3  μm, sample plane) used to perform cross-correlation measurements at higher resolution of the flow orientation (f) and of the speed (g). (h) CCFs calculated for each line of the grid [from top to bottom in (g)] between the circled spots and the corresponding downstream spots along the same line. The mean value of the flow speed retrieved from the slope of the time-of-flight as a function of the spot–spot distance is 1202±5  μm/s, with an accuracy of 0.4%. The acquisition frame rate in these experiments was equal to 986 fps. (i) Analysis of the dependence of the CCF amplitude as a function of the distance of the spots |R| normalized to the size of the spot ω0. The symbols refer to the A’ (filled squares), B’ (circles), C’ (up-triangles), and D’ (open squares) lines in (g).


The results shown in Fig. 5(b), respectively, obtained before (BC) and after (AC) the aberration correction, indicate that in both cases one of the CCFs (corresponding to the “C” spot) has dominant amplitude. The amplitude of the cross correlation was approximately threefold higher for the AC case (1.8) than for the BC case (0.6). We also found a systematically sharper correlation peak (60% reduction in width). We assign, therefore, the flow direction to the vector joining the corresponding couple of spots [from C’ to C in the example given in Fig. 5(a)].

The amplitude and the shape of the cross-correlation peak depend also on the threshold parameter, Th, and the size of the ROI over which the spot signal is averaged. The use of threshold values Th25 increases the amplitude of the cross-correlation peak of almost threefold [Fig. 5(c)]. As to the ROI size around each spot, the minimal choice m×m=1×1  pixels gives the largest amplitude of the CCF [Fig. 5(d)]. However, the result is also very sensitive to the exact positioning of the ROI on the irradiation spot. On the other hand, by increasing the ROI size above 5×5  pixels, the amplitude of the CCF decreases markedly [Fig. 5(d)]. Therefore, we set the threshold value at Th=25 and adopted a 3×3  pixels ROI in the following analyses. This choice will depend on the brightness of the dye labeling the tracers and the background of the flowing fluid.

We can monitor the flow speed in a wide range of values by measuring the time of flight from the maximum of the CCFs [Fig. 5(e)] between two spots. We tested it for three different values of the velocity (120, 400, and 1200  μm/s) retrieving the values set by the syringe actuator with an accuracy (the difference from the expected value) better than 3% and an uncertainty (variability of the measurement) of at most 6% (123±6, 395±25, and 1185±70  μm/s). The relative increase in the signal/noise obtained by correcting the aberrations is more and more effective the faster is the flow speed, as can be seen in Fig. 5(e) by comparing the dashed (corrected) to the solid (uncorrected) lines.

In order to refine the flow field mapping, we can device more elaborate illumination patterns with a larger number of spots. In this case, the correction of the aberrations was even more essential due to the rapid decrease of the two-photon excitation efficiency per spot with their number. An example on a 4×4 array is given in Fig. 5(f) (4×4 spots array, spacing Δx=7.0  μm, and Δy=5.3  μm in the sample plane). The analysis of the cross-correlation amplitudes (data not shown) indicates that the preferential flow direction is from the central spot to the “D” spot, within an accuracy angle determined by our choice of the array geometry that in the case of Fig. 5(f) θ5  μm21  μm12  deg. By means of a wider array, we can also reduce the uncertainty on the flow speed that can be measured for increasing distances between the spots [as sketched in Fig. 5(g)]. An example given in Fig. 5(h) suggests that by interpolating the time-of-flight time as a function of the spots distance the uncertainty on the speed can be reduced to at least 0.4% (speed=1202±5  μm/s).


Curved microchannel

The possibility offered by the SLM to create arbitrary patterns on the sample can be further exploited to characterize flow patterns. To illustrate this feature, we have employed a tightly curved microchannel [channel width 2R=700  μm; channel inner radius ρ=1225  μm; Fig. 6(a)] obtained by PDMS casting, sealed on a glass slide. The presence of uncured PDMS on the glass slide that was used to seal the microchannel, induced additional aberrations on the light wave front. We have computed the CCFs between couples of spots in a 5×2 array [Fig. 6(b), Δx=14.8  μm, Δy=4.6  μm, field of view (a) in Fig. 6(a)]. We found non-negligible cross-correlation amplitudes for all the couples of spots lying along the horizontal lines in Fig. 6(b). In this case, the correction of the aberrations is more effective than for the glass capillary: we could increase by a factor 4±0.3 the signal/background in the time-of-flight peak of the function by applying the aberrations-correction algorithm [Fig. 6(c)]. The flow speed uncertainty decreases almost threefold from 1.8% [v=900±17  μm/s, open squares, Fig. 6(c)] to 0.7% [v=920±7  μm/s, filled squares, Fig. 6(c)].

Fig. 6

(a) Fluorescence confocal image of horseshoe shaped microchannel employed to simulate a tightly curved capillary, filled with Rhodamine (the small white squares indicated with “a,” “b,” and “c” are the ROIs, in which the flow field was computed). (b) Maximum projection image of a 5×2 grid (interspot spacing 14.8×4.6  μm, sample plane) used to perform cross-correlation measurements of microbeads drag by a continuous flow. The field of view of this image is the white small open square “a” in (a). A similar image is obtained by the maximum projection of the stack of images acquired on the field of view “b” in (a). (c) CCFs computed between the couple of spots along line 3 in (b) computed on data obtained with (filled squares) and without (open squares) the correction of aberrations (frame rate=855  fps). The flow speed recovered from Gaussian fit was 920±7  μm/s and 900±17  μm/s   for these two cases, respectively. (d) CCFs computed between the pairs of spots laying on the lines 1 to 5 in the field of view “b” of (a), after the corrections of the aberrations. The speed retrieved from the Gaussian fit (solid lines) was 410±5  μm/s. (e) Maximum projection image of a semicircular disposition of the spots referring to the field of view “c” in (a). The cross-correlation functions between the central spot (H, in red circle) and each of the A–G spots on the circumference (orange circles) are reported in (f) (frame rate=986  fps). The speed retrieved by the best Gaussian fit (solid line) of the pair H–C of spots was 457±3  μm/s.


These results are confirmed also in different positions of the microchannel. For example by averaging over all the pairs of spots lying along the lines 1 to 5 in the field of view (b) [see Fig. 6(a)] that lies closer to the channel wall (at about 240  μm0.7R from the channel axis), we find v=408±5  μm/s [Fig. 6(d), solid lines], that agrees well with the estimate vmax[1(0.7RR)2]420  μm/s.

We have exploited the SLM also to shed a circular array of spots [shown as maximum projection image in Fig. 6(e)] on the sample plane at the position marked by the field of view (c) in Fig. 6(a). Such patterns could help in discriminating the direction of the flow for curved flow lines. The results of the cross-correlation measurements between the central spot and each of the nearest neighbor spots [Fig. 6(f)] indicate that the flow points toward the spot “C,” in agreement with the shape of the microchannel, with a best fit value of the speed 457±3  μm/s, not dissimilar from the speed in the field of view (b).




Aberrations Correction

Our algorithm is based on grids of spots shined and shifted regularly on the sample only by means of the SLM. The accuracy in the prediction of the spots position of a grid projected onto the sample plane is fundamental to achieve high-quality image reconstruction and accurate positioning of the observation volumes for the computation of the CCFs. We have shown that our algorithm, based on a combination of the BDA and image filtering, retrieves the spots position within 50  nm (well within the pixel size) for signal/background ratios as low as 1.35 [see Fig. 2(b)].

Our aberrations-correction algorithm is based on an intensity metric averaged over the spots of the grid [Eq. (6)] that spans the whole field of view. We can routinely reach a relative improvement of the signal-to-background ratio between the raw and the corrected spots between 3 [Fig. 2(c)] and 5 [Fig. 3(d)] over the whole field of view. However, the spot–spot spacing in the grid affects the improvement of the spot signal/background obtained by means of the aberrations correction. It decreases from 80±10 to 3±0.7 when the spot–spot distance decreases from 67 to 8  μm. Obviously, we would like to keep the spot–spot distance as small as possible, and the number N of spots in the array as large as possible in order to reduce the acquisition time. The improvement in the excitation efficiency obtained through the aberrations correction allows us to use arrays as large as 13×13 (see Fig. 4).

The corrections of the aberrations can be influenced also by (anti)correlation effects in the modes or in the image regions. Regarding the first possibility, we have tested our algorithm for the aberrations correction on a grid of spots for a possible cross-talk between pair of modes h and k. This parameter could be measured through the ratio of the out-of-trace terms, γk,hk of the γ matrix and a geometric average of k’th and the h’th modes term γk,kγh,h. However, we estimated experimentally the degree of coupling of the k’th mode with all the others, from the ratio of the metric obtained by correcting for the k’th mode only over the metric value obtained by correcting all the modes (GSD algorithm). In our experiments, this ratio was typically <0.9, indicating that the coupling between the modes was limited.

Regarding the second possibility, we should consider that the optical aberrations may differ substantially over the field of view. The contribution to the metric value of brighter spots is higher with respect to the dimmer ones that typically lie at the margin of the FOV: difference in the spots brightness can be serious for large (>100  μm)  FOV. Moreover, we expect that living tissues display anisoplanatic aberrations.47,48,55 We have tested the optical setup for possible anisoplanatism by performing corrections on arrays on adjacent ROIs covering a wide field of view (see Fig. 3). Operatively, on each ROI, we could compute the amplitudes of the Zernike polynomials that could be used to correct the optical aberrations locally and then reconstruct the wide field of view by stitching together the ROIs. However, a few general conclusions can also be drawn. The overall correction for each ROI shown in Fig. 3(a), measured as the average of all the Zernike amplitudes over the ROI, ζ=17n=511zn2, does not increase dramatically when moving from the center of the FOV (ζ0.6) to the periphery ζ0.7±0.2. Moreover, the most effective modes as measured by the mode amplitude averaged over the 9 ROIs, zn2ROI, are n=8, 5, and 10 (zn2ROI=0.7, 1.0, 1.2, respectively). This type of analysis can be in principle used to evaluate operatively the anisoplanatism over the desired field of view.

The interplay of the limited spatial frequency resolution and content of the SLM with the zero-order reflection produces unexpected effects in the aberration correction as we can judge from the analysis of the spot size and the signal/background ratio of the laser spots array reported in Figs. 3 and 4. In fact, we observed an apparent reduction of the spot size if measured on top of the background level [the full-width at half-maximum of the spots measured on top of the background decreases from 2.3±0.4  μm to 1.3±0.2  μm, Figs. 3(c)3(e)] when we decrease the spot–spot distance in the grid from 67 to 8  μm. This behavior, which is opposite to that observed for the signal/background level [see Fig. 3(c) inset], is an artifact related to the increase of the background level on the grids at low spot–spot spacing. This background increase is probably due to the discrete nature and the finite size56 of the hologram encoded on the SLM that prevents to accurately encode the hologram containing the Fresnel zone mask needed to correct for the zero-order reflection and a set of closely spaced spots.


Image Reconstruction and Aberrations Correction

Having tested the accuracy of the spot position retrieval and the aberration correction algorithm, we moved to combine them in the image reconstruction algorithm. The results shown in Fig. 4 indicate that modal corrections of the optical aberrations can be implemented in the SLM-based imaging algorithm based on the “shift and snap” method.31 The total acquisition time for the image shown in Fig. 4(e) (100×100  pixels, 400 images) is about 25 s (10-ms exposure time per frame).

The major limitation to the image acquisition rate comes from the refresh time of the SLM, 1/60  s for the Hamamatsu X10468-7 used here. The sample brightness would not limit the image acquisition rate substantially, unless for very dim samples. In fact, from a voxel350×350×1000  nm3 in a sample with a protein concentration 1  μM and a count per molecule 3000  counts/s (for a 10% collection efficiency), we should get 600  counts/ms or 240 electrons per pixel of an image acquired at 1 kHz (conversion factor 0.49  electrons/count for Hamamatsu Orcad camera), corresponding to an acceptable level of noise 6%. The minimum theoretical image acquisition time is then largely limited by the SLM refresh time and the synchronization between the SLM and the camera. In our case, this effective time is not <32  ms. For a 70×70  μm2 field of view, 0.14  μmpixel size, 512×512  pixels (corresponding to about 105 frames) the minimum acquisition time is then 9 min.


Aberrations Correction and Cross-Correlation Spectroscopy

The increase of the efficiency of the nonlinear excitation obtained by the optical aberrations correction, the accuracy in the laser spot positioning, and the use of no movable parts make our setup particularly suited to characterize flows in complex vessel networks,57,58 in confined environments,59 or in vivo8 by means of cross-correlation methods. For these applications, it is essential to be able to select the observation volumes on an image collected with the same setting of the optical setup, particularly for low brightness and signal/noise ratio samples.31 This can be easily achieved by means of our SLM-based microscope as we verify here by studying the flow of microbeads (1-μm size) at decreasing values of the excitation intensity as reported in Fig. 5 for the case of a straight glass microchannel. The amplitude of the cross correlation was 3 to 4 times higher for the aberration corrected CCFs than for the raw data functions [see Figs. 5(b) and 6(c)].

Our spots-array method can be tailored easily to a variety of shapes of the vessels and vessel networks as shown in principle here for a curved microchannel (Fig. 6). The uncertainty in the estimate of the flow direction is determined by the choice of the grid. In particular, for a rectangular grid, the relevant parameter is the angular size with which the upstream spot sees the downstream spot. In Fig. 5(a), this angle (for example between the A and A’ spots) is θ±1  μm4.6  μm±12  deg. The resolution is instead determined by the spacing of the spots along a direction perpendicular to the flow [for example between A and B in Fig. 5(a)], that is 4.6  μm14.8  μm17  deg. Given the results in Figs. 3(c)3(e), the interspot distance along a vertical line in Fig. 5(a) cannot be reduced much without decreasing substantially the signal/background ratio. However, one could increase the distance between the vertical lines of spots or shift vertically one column of spots with respect to the other: in both cases, we can easily halve the angular resolution with respect to Fig. 5(a).

Moreover, we can gain additional information on the tracer and on the flow from the decrease of the amplitude of the CCF as a function of the distance |R| between the spots. As it can be seen in Fig. 5(h), the CCF amplitude stays either approximately constant or decreases slightly with this distance. This trend can be due to a low Peclet number60 (i.e., large diffusion coefficients D such that |v||R|<D) that would imply an increase in the width of the CCF peak and a parallel decrease of its amplitude. Alternatively, the amplitude of the CCF peak decreases when the flow line is not perfectly aligned with the vector R, and therefore, α0 [Eq. (14)]. A simple analysis of Eq. (14) suggests that the amplitude decay should follow an exponential law exp(α2/β2), where β=ω0/|R| is the angle under which the upstream spot sees the downstream spot (whose size is ω0). An analysis of the trend of the ln[GX12(τmax)] as a function of 1β2=|R|2/ω02 provides misalignment angles at most of α5  deg [for the A’ line in Fig. 5(g)], well within the uncertainty determined by the array geometry. In the cases investigated here, no systematic correlation of the width of the CCF peak with the estimated α value was found, indicating that the Peclet number in our case is large (in fact we estimate Pc=|v||R|/D19,000).

We notice also that the possibility to create arbitrary patterns of laser spots in the sample would allow us to accurately test the flow field pattern. In the case of the curved microchannel reported in Fig. 6(a) (width 2R=700  μm and inner radius ρ=1225  μm), we could have deviations from the axial flow such as Dean flow60 that could be detected in the image plane. However, since the Reynold number in our experiments is relatively low, Ry0.3, any secondary Dean flow60 in the radial direction is within the current uncertainty on the flow velocity direction. In fact, we expect45 a radial component of the velocity vrRyRρ|v|, or an angle of misalignment δθRyRρ with respect to the laminar flow. In our case, this value δθ4.9  deg, well within the uncertainty due to the grid spacing, θ±12  deg. However, as discussed above, this limitation can be overcome by means of an appropriate choice of the spots array.

The time resolution of the flow speed measurements presented here depends on two factors: the concentration of the tracer objects and the uncertainty that we want to reach on the flow speed. We typically acquired traces of 3000 to 4000 frames at a frame rate of 986 Hz for a total acquisition time of about 3 or 4 s. With the particle density of this study, the amplitude of the cross-correlation peak is fairly constant with an error smaller than 10% for total acquisition time 1  s. The spread in the CCF amplitude raises to 50% for a total acquisition time 0.4  s. However, the lagtime position of the peak has an uncertainty smaller than 3% for total duration of the order of 0.4 s. Therefore, we can measure the flow speed within a few percent points with a time resolution of 0.4 s. These considerations could vary with microbeads of different brightness and different concentrations.

Finally, we notice that the correction of the aberrations positively affects the uncertainty of the flow speed measurement. For the case of a flow in a straight microchannel [Fig. 5(b)], the width of the CCF peak decreases almost twofold when correcting the aberrations (from δτ(0)=2.5±0.2 to δτ=1.6±0.1  ms). The corresponding lagtime does not change appreciably with the application of the corrections but its uncertainty decreases, passing from τmax(0)=9.7±0.09  ms to τmax=9.6±0.03  ms. From these values, we can estimate the flow speed (the interspot distance is d=14.7±0.3  μm) that results v(0)=dτmax(0)=1526±14  μm/s, when no aberrations correction is applied, and v=1542±6  μm/s when the aberrations are corrected, with a twofold reduction of the uncertainty from 0.8% to 0.3%. The width of the CCF peak provides us with an indirect measurement of the beam spot radius ω0 through ω0vδτ [see Eq. (14)]. The values of the spot size for the two cases are ω04.0±0.3  μm and 2.5±0.1  μm, for the two cases when we do not correct for the aberrations and when we correct for them, respectively. These figures agree with the measurement reported on fixed grids in Figs. 4(a) and 4(b). The twofold decrease in the speed uncertainty, obtained when applying the aberrations correction, is found also at various speed values [see Fig. 5(e)]: the uncertainty reduces of 0.5±0.2 times in the range 120v1200  μm/s. Moreover, the accuracy increases threefold for the absolute speed value: for a nominal flow speed v=120  μm/s, we find v=121.4±0.06  μm/s (accuracy 1%) when the aberrations are corrected, and v(0)=124.0±0.1  μm/s (accuracy 3%) when no aberrations correction is applied [Fig. 5(e)].



By means of a purely diffractive microscope based on an SLM, we collected high resolution, aberration corrected two-photon images, on which we positioned with high-accuracy arrays of laser spots and compute CCFs between couples of them. To this purpose, we have combined the shift-and-snap algorithm based on the synchronization of the SLM and the camera with AO minimization of an intensity-based metric. The inclusion of the SLM in the optical path allows also to shed arbitrary light patterns on the sample and to perform cross-correlation measurements over multiple pairs of spots, therefore, expanding the possibility of previous setups.61

In our experience, even for controlled samples in vitro, such as microfluidic devices obtained by PDMS casting, the correction of the aberrations is essential to reach good contrast in the CCFs and to measure accurately the flow velocity. We obtained with the present setup a routine fivefold improvement of the signal/background on the image and an increase in the contrast of the CCF that can reach three times [Figs. 5(b) and 6(c)]. This improvement result in a twofold decrease of the uncertainty in the evaluation of the flow speed and, at least for flow speeds <1000  μms, a threefold increase in the accuracy of the flow speed evaluation.

We expect that the features displayed by the SLM-based correlation microscope can be valuable for the study of hemodynamics in small animal models, such as Zebrafish embryos. In this case, even for transgenic albino embryonic lines, the scattering from the tissue will induce aberrations, possibly anisoplanatic, of the wave front. The possibility of our setup to operate separate aberrations correction over different regions of interest in the field of view will be a valuable help in collecting high-quality dynamic data.


The authors have no relevant financial interests in this article and no potential conflicts of interest to disclose.


A small part of this work has been previously published in the conference proceeding “Scanless nonlinear optical microscope for image reconstruction and space-time correlation analysis,” N. Ceffa; F. Radaelli; P. Pozzi; M. Collini; L. Sironi; L. D’alfonso; G. Chirico, Proc. SPIE. 10333, Optical Methods for Inspection, Characterization, and Imaging of Biomaterials III. We acknowledge the support of the funding 2015-ATE-0007 by the University of Milano-Bicocca to G. C. We are grateful to Dr. Paolo Pozzi (TU-Delft, NL) for helpful discussions about the AO algorithms and to Prof. F. Auricchio and Dr. S. Marconi of the university of Pavia, Department of Civil Engineering and Architecture for the 3D printing of the molds.



C. Cierpka and C. J. Kaehler, “Particle imaging techniques for volumetric three-component (3D3C) velocity measurements in microfluidics,” J. Visual., 15 (1), 1 –31 (2012). Google Scholar


H. Kim, J. Westerweel and G. E. Elsinga, “Advances in 3D velocimetry,” Meas. Sci. Technol., 24 (2), 024007 (2013). MSTCEP 0957-0233 Google Scholar


M. Bouzin et al., “An intermittent model for intracellular motions of gold nanostars by k-space scattering image correlation,” Biophys. J., 109 (11), 2246 –2258 (2015). BIOJAU 0006-3495 Google Scholar


F. Helmchen and W. Denk, “Deep tissue two-photon microscopy,” Nat. Methods, 2 (12), 932 –940 (2005). 1548-7091 Google Scholar


M. Caccia et al., “Image filtering for two-photon deep imaging of lymphonodes,” Eur. Biophys. J. Biophys. Lett., 37 (6), 979 –987 (2008). Google Scholar


M. J. Booth, “Adaptive optics in microscopy,” Optical and Digital Image Processing: Fundamentals and Applications, 295 –322 Wiley-VCH Verlag GmbH & Co., Weinheim (2011). Google Scholar


S. Gharbi et al., “Reduction of chromatic dispersion using multiple carrier frequency patterns in SLM-based microscopy,” Appl. Opt., 56 (23), 6688 –6693 (2017). APOPAI 0003-6935 Google Scholar


M. Hasler et al., “Modulator-based phase contrast microscopy,” Opt. Eng., 54 (4), 043107 (2015). Google Scholar


C. Lingel, T. Haist and W. Osten, “Spatial-light-modulator-based adaptive optical system for the use of multiple phase retrieval methods,” Appl. Opt., 55 (36), 10329 –10334 (2016). APOPAI 0003-6935 Google Scholar


N. Matsumoto et al., “Aberration correction considering curved sample surface shape for non-contact two-photon excitation microscopy with spatial light modulator,” Sci. Rep., 8 9252 (2018). SRCEC3 2045-2322 Google Scholar


M. A. A. Neil et al., “Adaptive aberration correction in a two-photon microscope,” J. Microsc., 200 105 –108 (2000). JMICAR 0022-2720 Google Scholar


D. Débarre et al., “Image-based adaptive optics for two-photon microscopy,” Opt. Lett., 34 (16), 2495 –2497 (2009). OPLEDP 0146-9592 Google Scholar


C. Wang et al., “Multiplexed aberration measurement for deep tissue imaging in vivo,” Nat. Methods, 11 1037 –1040 (2014). 1548-7091 Google Scholar


J. Mertz, H. Paudel and T. G. Bifano, “Field of view advantage of conjugate adaptive optics in microscopy applications,” Appl. Opt., 54 (11), 3498 –3506 (2015). APOPAI 0003-6935 Google Scholar


O. Albert et al., “Smart microscope: an adaptive optics learning system for aberration correction in multiphoton confocal microscopy,” Opt. Lett., 25 (1), 52 –54 (2000). OPLEDP 0146-9592 Google Scholar


S. Bonora et al., “Wavefront correction and high-resolution in vivo OCT imaging with an objective integrated multi-actuator adaptive lens,” Opt. Exp., 24 (17), 21931 –21941 (2015). OPEXFF 1094-4087 Google Scholar


F. Gao et al., “3D Imaging of flow patterns in an internally-pumped microfluidic device: redox magnetohydrodynamics and electrochemically-generated density gradients,” Anal. Chem., 85 (9), 4414 –4422 (2013). ANCHAM 0003-2700 Google Scholar


L. Sironi et al., “In vivo flow mapping in complex vessel networks by single image correlation,” Sci. Rep., 4 7341 (2014). SRCEC3 2045-2322 Google Scholar


W. Krieger et al., “Imaging fluorescence (cross-) correlation spectroscopy in live cells and organisms,” Nat. Protoc., 10 (12), 1948 –1974 (2015). 1754-2189 Google Scholar


N. G. Ceffa et al., “Spatiotemporal image correlation analysis for 3D flow field mapping in microfluidic devices,” Anal. Chem., 90 (1), 2277 –2284 (2018). ANCHAM 0003-2700 Google Scholar


C. A. Marquezin et al., “Image cross-correlation analysis of time varying flows,” Anal. Chem., 88 (14), 7115 –7122 (2016). ANCHAM 0003-2700 Google Scholar


K. Brinkmeier, J. S. Doerre and M. Eigen, “Two-beam cross-correlation: a method to characterize transport phenomena in micrometer-sized structures,” Anal. Chem., 71 (3), 609 –616 (1999). ANCHAM 0003-2700 Google Scholar


P. S. Dittrich and P. Schwille, “Spatial two-photon fluorescence crosscorrelation spectroscopy for controlling molecular transport in microfluidic structures,” Anal. Chem., 74 (17), 4472 –4479 (2002). ANCHAM 0003-2700 Google Scholar


T. Dertinger et al., “Two-focus fluorescence correlation spectroscopy: a new tool for accurate and absolute diffusion measurements,” Chem. Phys. Chem., 8 (3), 433 –443 (2007). Google Scholar


C. B. Mueller et al., “Precise measurement of diffusion by multi-color dual-focus fluorescence correlation spectroscopy,” Eur. Phys. Lett., 83 (4), p1 –p5 (2008). Google Scholar


M. J. Rossow et al., “Raster image correlation spectroscopy in live cells,” Nat. Protoc., 5 (11), 1761 –1774 (2010). 1754-2189 Google Scholar


L. Sacconi et al., “Optical recording of electrical activity in intact neuronal net-works with random access second-harmonic generation microscopy,” Opt. Exp., 16 14910 –14921 (2008). OPEXFF 1094-4087 Google Scholar


P. A. Kirkby, K. M. Srinivas Nadella and R. A. Silver, “A compact acousto-optic lens for 2D and 3D femtosecond based 2-photon microscopy,” Opt. Exp., 18 13721 –13745 (2010). OPEXFF 1094-4087 Google Scholar


R. Salomé et al., “Ultrafast random-access scanning in two-photon microscopy using acousto-optic deflectors,” J. Neurosci. Methods, 154 161 –174 (2006). JNMEDT 0165-0270 Google Scholar


V. Nikolenko et al., “SLM microscopy: scanless two-photon imaging and photo-stimulation with spatial light modulators,” Front. Neural Circuits, 2 5 (2008). Google Scholar


P. Pozzi et al., “High-throughput spatial light modulation two-photon microscopy for fast functional imaging,” Neurophotonics, 2 (1), 015005 (2015). Google Scholar


D. Gandolfi et al., “The spatio temporal organization of cerebellar network activity resolved by two-photon imaging of multiple single neurons,” Front. Cell. Neurosci., 8 92 (2014). Google Scholar


R. W. Gerchberg and W. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik, 35 (2), 237 –246 (1972). OTIKAJ 0030-4026 Google Scholar


R. Di Leonardo, F. Ianni and G. Ruocco, “Computer generation of optimal holograms for optical trapping,” Opt. Express, 15 (4), 1913 –1922 (2007). OPEXFF 1094-4087 Google Scholar


H. Zhang et al., “Elimination of a zero-order beam induced by a pixelated spatial light modulator for holographic projection,” Appl. Opt., 48 (30), 5834 –5841 (2009). APOPAI 0003-6935 Google Scholar


J. M. Girkin, S. Poland and A. J. Wright, “Adaptive optics for deeper imaging of biological samples,” Curr. Opin. Biotechnol., 20 (1), 106 –110 (2009). CUOBE3 0958-1669 Google Scholar


M. J. Booth, “Adaptive optical microscopy: the ongoing quest for a perfect image,” Light Sci. Appl., 3 e165 (2014). Google Scholar


W. Moon et al., “Computer-aided tumor detection based on multi-scale blob detection algorithm in automated breast ultrasound images,” IEEE Trans. Med. Imaging, 32 1191 –1200 (2013). ITMID4 0278-0062 Google Scholar


V.-G. Yalla, W. Su and L. G. Hassebrook, “Multi-spot projection, tracking and calibration,” Proc. SPIE, 5106 221 –232 (2003). PSISDG 0277-786X Google Scholar


A. Jesacher, M. J. Booth, “Sensorless adaptive optics for microscopy,” Adaptive Optics for Biological Imaging, 177 –189 CRC Press, Boca Raton, Florida (2013). Google Scholar


M. J. Booth, “Wavefront sensorless adaptive optics for large aberrations,” Opt. Lett., 32 (1), 5 –7 (2007). OPLEDP 0146-9592 Google Scholar


M. Born and E. Wolf, Principles of Optics, 2nd ed.Pergamon Press, Bath, UK (1964). Google Scholar


V. N. Mahajan and G.-M. Dai, “Orthonormal polynomials in wavefront analysis: analytical solution,” J. Opt. Soc. Am. A, 24 (9), 2994 –3016 (2007). JOAOD6 0740-3232 Google Scholar


R. J. Noll, “Zernike polynomials and atmospheric turbulence,” J. Opt. Soc. Am., 66 (3), 207 –211 (1976). JOSAAH 0030-3941 Google Scholar


A. Thayil and M. J. Booth, “Self calibration of sensorless adaptive optical microscopes,” J. Eur. Opt. Soc., 6 11045 (2011). 1990-2573 Google Scholar


J. Zeng et al., “3D resolved mapping of optical aberrations in thick tissues,” Biomed. Opt. Exp., 3 (8), 1898 –1913 (2012). BOEICL 2156-7085 Google Scholar


G. Han et al., “Isoplanatic patch of the human eye for arbitrary wavelengths,” Opt. Commun., 410 (1), 811 –816 (2018). OPCOB8 0030-4018 Google Scholar


R. Aviles-Espinosa et al., “Measurement and correction of in vivo sample aberrations employing a nonlinear guide-star in two-photon excited fluorescence microscopy,” Biomed. Opt. Express, 2 (11), 3135 –3149 (2011). BOEICL 2156-7085 Google Scholar


D. Débarre et al., “Adaptive optics for structured illumination microscopy,” Opt. Exp., 16 (13), 9290 –9305 (2008). OPEXFF 1094-4087 Google Scholar


W. H. Press et al., “Parabolic interpolation and Brent’s method in one dimension,” Numerical Recipes: The Art of Scientific Computing, Cambridge University Press(2007). Google Scholar


L. V. S. Boyd, Convex Optimization, 463 –475 Cambridge University Press, New York (2004). Google Scholar


M. A. Digman and E. Gratton, “Lessons in fluctuations correlation spectroscopy,” Annu. Rev. Phys. Chem., 62 645 –668 (2011). ARPLAP 0066-426X Google Scholar


S. Chen, P. Tartaglia and P. Pusey, “Light-scattering from independent particles—non-Gaussian correction to clipped intensity correlation-function,” J. Phys. A Math. Gen., 6 (4), 490 –495 (1973). Google Scholar


K. Bacia, S. A. Kim and P. Schwille, “Fluorescence cross-correlation spectroscopy in living cells,” Nat. Methods, 3 (2), 83 –89 (2006). 1548-7091 Google Scholar


M. Laslandes et al., “Increasing the field of view of adaptive optics scanning laser ophthalmoscopy,” Biomed. Opt. Exp., 8 (11), 4811 –4826 (2017). BOEICL 2156-7085 Google Scholar


M. K. Kim, “Principles and techniques of digital holographic microscopy,” SPIE Rev., 1 (1), 018005 (2010). Google Scholar


N. G. Ceffa et al., “Spatiotemporal image correlation analysis of blood flow in branched vessel networks of Zebrafish embryos,” J. Biomed. Opt., 22 (10), 106008 (2017). JBOPFO 1083-3668 Google Scholar


L. Fieramonti et al., “Quantitative measurement of blood velocity in Zebrafish with optical vector field tomography,” J. Biophotonics, 8 (1-2), 52 –59 (2015). Google Scholar


D. J. Rowland, H. H. Tuson and J. S. Biteen, “Resolving fast, confined diffusion in bacteria with image correlation spectroscopy,” Biophys. J., 110 (10), 2241 –2251 (2016). BIOJAU 0006-3495 Google Scholar


T. M. Squires and S. R. Quake, “Microfluidics: fluid physics at the nanoliter scale,” Rev. Mod. Phys., 77 (3), 977 –1026 (2005). RMPHAT 0034-6861 Google Scholar


P. Pozzi et al., “Electron multiplying charge-coupled device-based fluorescence crosscorrelation spectroscopy for blood velocimetry on zebrafish embryos,” J. Biomed. Opt., 19 (6), 067007 (2014). JBOPFO 1083-3668 Google Scholar

Biographies of the authors are not available.

© The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Maddalena Collini, Fabrizio Radaelli, Laura Sironi, Nicolo G. Ceffa, Laura D'Alfonso, Margaux Bouzin, and Giuseppe Chirico "Adaptive optics microspectrometer for cross-correlation measurement of microfluidic flows," Journal of Biomedical Optics 24(2), 025004 (27 February 2019).
Received: 11 September 2018; Accepted: 4 December 2018; Published: 27 February 2019

Back to Top