*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.

## 1.

## Introduction

The measurement of flows in microfluidic systems or in blood vessels can be obtained either by means of single particle tracking^{1}^{–}^{3} 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 tissues^{4} 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 state^{4}^{,}^{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 modality^{7} 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 correction^{13}^{,}^{14}. In the first case, the wave front is described as the superposition of $\mathrm{\Xi}$ 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 wavelength^{15} 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 analysis^{17} to flow mapping on single images.^{18} It has recently acquired great impact in cellular biology, where it is applied on 2-D images^{19}^{,}^{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 sample^{22}^{–}^{25} 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 deflector^{27}^{–}^{29} or an SLM^{30} 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 neurons^{32} 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.

## 2.

## Materials and Methods

## 2.1.

### 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 $\lambda =800\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$. The polarization of the laser light is perpendicular to the optical table. The laser is collimated and expanded by a $3.5\times $ beam expander (${f}_{1}=40\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ and ${f}_{2}=140\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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\pi $ on the beam wave front with $792\times 600\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}$ resolution ($20\text{-}\mu \mathrm{m}\text{\hspace{0.17em}}\text{pixel}$ size) and 8-bit pixel depth. The SLM plane is conjugated through a $0.65\times $ beam reducer (${f}_{3}=200\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ and ${f}_{4}=125\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$) to the entrance pupil of a microscope objective ($20\times $ 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}

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 (${f}_{\mathrm{LT}}=150\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$) on a Hamamatsu CMOS Orca Flash 4.0 Camera (pixel size $6.5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$).

## 2.2.

### 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.

## 2.3.

### Grid Projection

A grid of $N\times 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\pi $ 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\times 600\text{\hspace{0.17em}\hspace{0.17em}}\text{pixel}$). We calculate two conversion factors ${\kappa}_{\mathrm{col}}$ and ${\kappa}_{\text{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)].

## 2.4.

### 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\epsilon |$, $\epsilon <0$, and we add a quadratic phase pattern (FL) to the GSW-phase.^{31} The FL has the following expression:

## Eq. (1)

$${\varphi}_{\mathrm{FL}}(\overrightarrow{x})=2\pi \kappa \frac{\epsilon}{{f}^{2}}{|\overrightarrow{x}|}^{2},$$## 2.5.

### 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.

## 2.5.1.

#### 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.

## 2.5.2.

#### 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 model^{39} 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 $({x}_{s},{y}_{s})$ and the camera domain $({x}_{c},{y}_{c})$ can be expressed as

## Eq. (2)

$${x}_{c}={({a}_{xx}{x}_{s}+{b}_{xx})}^{2}{({a}_{xy}{y}_{s}+{b}_{xy})}^{2},\phantom{\rule{0ex}{0ex}}{y}_{c}={({a}_{yx}{x}_{s}+{b}_{yx})}^{2}{({a}_{yy}{y}_{s}+{b}_{yy})}^{2}.$$We adopt a matrix notation, where $\mathbf{P}$ is a $N\times 2$ matrix ($N$ is the number of spots in the grid) whose rows contains the coordinates $({x}_{c,i},{y}_{c,i})$ of the $i$’th spot in the camera domain, $\mathbf{C}$ is a $N\times 9$ matrix whose $i$’th row contains the products ${\{{{x}_{s,i}}^{n}{{y}_{s,i}}^{m}\}}_{m,n=\mathrm{0,1},2}$, and $\mathbf{A}$ is a $9\times 2$ matrix that contains, in each row, the expansion coefficients of ${x}_{c,i}$ and ${y}_{c,i}$ into polynomials [Eq. (2)]. Once the matrices $\mathbf{P}$ and $\mathbf{C}$ are obtained from the mapping calibration procedure, the transformation matrix $\mathbf{A}$ is simply calculated as^{39}

## Eq. (3)

$$\mathbf{A}={({\mathbf{C}}^{\mathrm{T}}\mathbf{C})}^{-1}{\mathbf{C}}^{\mathrm{T}}\mathbf{P}.$$The matrix $\mathbf{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 $\mathbf{A}$ matrix, the position of the new spots in the camera domain could be estimated as

## 2.6.

### 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.

## 2.6.1.

#### 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:

where ${\psi}_{i}$ is the $i$’th basis function and ${z}_{i}$ 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 modes

^{43}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 studies^{11}^{,}^{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 $\hat{{S}_{n}}$ integrated over ROIs (size $m\times m$) centered on each of the $N$ spots of the grid (${\{{x}_{c,i},\text{\hspace{0.17em}\hspace{0.17em}}{y}_{c,i}\}}_{i=1,\dots ,N}$) with the background measured by calculating the mean integrated intensity $\hat{B}$ over four ROIs (size $m\times 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

## 2.6.2.

#### Metric optimization

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

## Eq. (7)

$$M\approx q(1-\sum _{k=1}^{\mathrm{\Xi}}\sum _{l=1}^{\mathrm{\Xi}}{\gamma}_{kl}{z}_{k}{z}_{l})=(1-{\mathbf{z}}^{\mathrm{T}}\mathbf{\gamma}\mathbf{z}).$$In Eq. (7), $q$ and $\mathbf{\gamma}={\{{\gamma}_{kl}\}}_{k,\mathrm{l}=1\dots S}$ 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 $\mathbf{\gamma}$ is not strictly diagonal, the optimization of the metric is not straightforward. It is possible to diagonalize the matrix^{17}^{,}^{49} $\mathbf{\gamma}$ 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 (${z}_{l}=0$ if $k\ne l$). 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 algorithm^{51} 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.

## 2.7.

### 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 ${N}_{\text{row}}\times {N}_{\mathrm{col}}$ 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)

$${\varphi}_{P}(\overrightarrow{x})=2\pi {\overrightarrow{\delta}}_{P}\xb7\overrightarrow{x},$$The number of shift movements along horizontal and vertical directions (respectively, ${\mathrm{\Lambda}}_{H}$ and ${\mathrm{\Lambda}}_{V}$) is determined by the number of spots and the size in pixels of the reconstructed image:

## Eq. (9)

$${\mathrm{\Lambda}}_{H}=\mathrm{int}(\frac{W}{{N}_{\text{col}}}-1)\phantom{\rule[-0.0ex]{1em}{0.0ex}}{\mathrm{\Lambda}}_{V}=\mathrm{int}(\frac{H}{{N}_{\text{row}}}-1),$$## Eq. (10)

$${\mathrm{\Delta}}_{\mathrm{col}}^{\mathrm{SLM}}=\text{round}\left\{\frac{\mathcal{M}\xb7{\delta}_{\mathrm{col}}({\mathrm{\Lambda}}_{H}+1)}{{P}_{x}\xb7{\kappa}_{\mathrm{col}}}\right\}\phantom{\rule[-0.0ex]{2em}{0.0ex}}{\mathrm{\Delta}}_{\text{row}}^{\mathrm{SLM}}=\text{round}\left\{\frac{\mathcal{M}\xb7{\delta}_{\text{row}}({\mathrm{\Lambda}}_{V}+1)}{{P}_{y}\xb7{\kappa}_{\text{row}}}\right\},$$## Eq. (11)

$${\delta}_{\mathrm{col}}^{\mathrm{SLM}}=\frac{{\mathrm{\Delta}}_{\mathrm{col}}^{\mathrm{SLM}}}{{\mathrm{\Lambda}}_{H}+1}\phantom{\rule[-0.0ex]{2em}{0.0ex}}{\delta}_{\text{row}}^{\mathrm{SLM}}=\frac{{\mathrm{\Delta}}_{\text{row}}^{\mathrm{SLM}}}{{\mathrm{\Lambda}}_{V}+1}.$$Finally, the signal read from the spot in the $i$’th column (from left to right, $i=0,1,\dots ,{N}_{\mathrm{col}}$) and $j$’th row (from top to bottom, $j=\mathrm{0,1},\dots ,{N}_{\text{row}}$) of a grid that has been shifted $l$ times along the horizontal direction ($l=\mathrm{0,1},\dots ,{\mathrm{\Lambda}}_{H}$) and $k$ times along the vertical direction ($k=\mathrm{0,1},\dots ,{\mathrm{\Lambda}}_{V}$) is assigned in the synthetic reconstructed image to the pixel with indices $\{I,J\}$:

When we reconstruct an image, we usually integrate the signal from a $3\times 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 $({\mathrm{\Delta}}_{\mathrm{col}}^{\mathrm{SLM}},{\mathrm{\Delta}}_{\text{row}}^{\mathrm{SLM}})$ and the shift $({\delta}_{\mathrm{col}}^{\mathrm{SLM}},{\delta}_{\text{row}}^{\mathrm{SLM}})$ parameters can be converted in the sample plane by means of equivalent parameters called $(\mathrm{\Delta}x,\mathrm{\Delta}y)$ and $(\delta x,\delta y)$ used in Sec. 3.

## 2.8.

### 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)

$${G}_{X12}(\tau )=\frac{{\u27e8\delta {I}_{1}(t)\xb7\delta {I}_{2}(t+\tau )\u27e9}_{t}}{{\u27e8{I}_{1}(t)\u27e9}_{t}\xb7{\u27e8{I}_{2}(t)\u27e9}_{t}},$$^{4}

^{,}

^{6}

^{,}

^{52}

## Eq. (14)

$${G}_{X12}(\tau ;{\tau}_{D},{\tau}_{v})={G}_{\mathrm{Diff}\text{\hspace{0.17em}\hspace{0.17em}}}(\tau ;{\tau}_{D})\xb7\mathrm{exp}\{-\frac{{|\overrightarrow{v}|}^{2}}{{\omega}_{0}^{2}}\frac{[{\tau}^{2}+{\tau}_{v}^{2}-2\tau {\tau}_{v}\text{\hspace{0.17em}}\mathrm{cos}(\alpha )]}{1+\tau /{\tau}_{D}}\},\phantom{\rule{0ex}{0ex}}{G}_{\mathrm{Diff}}(\tau ;{\tau}_{D})=\frac{1}{\u27e8N\u27e9}\frac{1}{\sqrt{1+\tau /{\tau}_{D}}\sqrt{1+\xi \tau /{\tau}_{D}}}.$$In Eq. (14), ${\tau}_{D}={\omega}_{0}^{2}/4D$ is the diffusion time, $D$ is the diffusion coefficient of the colloids, ${\omega}_{0}$ is the beam waist, $\xi \cong 10-20$, accounts for the volume elongation along the optical axis, and ${\tau}_{v}=|\overrightarrow{R}|/|\overrightarrow{V}|$ is the time of flight between the observation volumes. When the flow axis is collinear ($\alpha =0$) to the interspot distance $\overrightarrow{R}$, Eq. (14) can be approximated by a simple Gaussian function in the lag time with maximum at the lag time ${\tau}_{\mathrm{max}}={\tau}_{v}$, for small values of the times ratio $\mathrm{\Gamma}={\tau}_{v}/{\tau}_{D}$. Numerical simulations proved^{6} that for $\mathrm{\Gamma}<0.1$ the discrepancy between ${\tau}_{v}$ and ${\tau}_{\mathrm{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 ${\tau}_{v}$ as ${\tau}_{\mathrm{max}}\cong {\tau}_{v}\text{\hspace{0.17em}}\mathrm{cos}(\alpha )$. Moreover, the amplitude ${G}_{X12}({\tau}_{\mathrm{max}};{\tau}_{D},{\tau}_{v})$ of the CCF decreases as a function of the noncollinearity angle $\alpha $ as ${G}_{X12}({\tau}_{\mathrm{max}};{\tau}_{D},{\tau}_{v})\approx \mathrm{exp}(-{\alpha}^{2}/{\beta}^{2})$, where $\beta $ 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 $\alpha $.

Equation (14) provides only the flow speed, through the measure of ${\tau}_{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 ${\tau}_{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.

## 3.

## Results

## 3.1.

### 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\times 5$ grid (spot–spot spacing $21.1\times 15.8\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ in the sample plane). In this example, the signal-to-background ratio was $9.5\pm 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 $\text{error}\cong \text{\hspace{0.17em}\hspace{0.17em}}50\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}\cong \frac{1}{7}\text{\hspace{0.17em}\hspace{0.17em}}\text{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.

## 3.2.

### Aberrations Correction

## 3.2.1.

#### 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 $\ge 20\%$. The improvement in the two-photon image of a $5\times 5$ grid (interspots spacing $28.8\times 21.8\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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 ($\text{signal}/\text{background}=7.3\text{\hspace{0.17em}}\pm 0.6$) and the corrected spots ($\text{signal}/\text{background}=2.4\pm 0.1$) lying on the dashed line in Fig. 2(c) is $3.0\pm 0.3$.

## 3.2.2.

#### Anisoplanatism

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\times 5$ grid in different regions of the FOV onto the Rhodamine slab [Fig. 3(a)]. The correction amplitude of each mode [${z}_{n}$ in Eq. (5), $5\le n\le 11$] 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.

## 3.2.3.

#### 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\pm 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\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, the signal/noise ratio on the corrected image decreases from $80\pm 10$ to $3\pm 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\pm 0.8\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ to $2.67\pm 0.02\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ for the same number of spots in the matrix [$5\times 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 $\text{distance}=27.8\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$. By changing the grid from a $3\times 3$ to a $7\times 7$ matrix, the full-width at half-maximum decreases from $2.90\pm 0.03$ to $2.33\pm 0.02\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ [Fig. 4(a)].

## 3.3.

### 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 $\lambda =800\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$) of the sample. Figure 4(c) shows an example of the images obtained before the aberrations correction, by means of a $13\times 13$ spots grid (spacing $\mathrm{\Delta}x=\mathrm{\Delta}y=27.4\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ in the sample plane) shifted with steps $\delta x=\delta y=702\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$ in the sample plane. The grid shift corresponds to $\sim 2$ CCD pixels in the camera plane. The grid rectangular distortion was quantified and corrected with correction factors ${\kappa}_{\mathrm{col}}=3.0\pm 0.05$ and ${\kappa}_{\text{row}}=4.0\pm 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\times 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\times 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\pm 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\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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\pm 0.1\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$. High-resolution imaging of the same sample features with a Leica SP5 confocal microscope (resolution $330\pm 30\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$) gives a width for these structures of $1.4\pm 0.5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$. These values indicate an optical resolution of the present setup of $1.5\pm 0.5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$.

## 3.4.

### 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\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}\times 800\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ size) and a horseshoe-shaped microchannel ($500\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}\times 500\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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\times 3\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$ area). The principle of cross-correlation spectroscopy^{52}^{,}^{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.

## 3.4.1.

#### 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\times N$ grid, we reduce the excitation intensity on each spot of the grid by a factor scaling as ${N}^{4}$. So, as a first choice, a $5\times 2$ rectangular grid of spots ($\mathrm{\Delta}x=14.8\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ and $\mathrm{\Delta}y=4.6\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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)].

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 ($\cong 1.8$) than for the BC case ($\cong 0.6$). We also found a systematically sharper correlation peak ($\cong 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 $\mathrm{Th}\ge 25$ 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\times m=1\times 1\text{\hspace{0.17em}\hspace{0.17em}}\text{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\times 5\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$, the amplitude of the CCF decreases markedly [Fig. 5(d)]. Therefore, we set the threshold value at $\mathrm{Th}=25$ and adopted a $3\times 3\text{\hspace{0.17em}\hspace{0.17em}}\text{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\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{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\pm 6$, $395\pm 25$, and $1185\pm 70\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{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\times 4$ array is given in Fig. 5(f) ($4\times 4$ spots array, spacing $\mathrm{\Delta}x=7.0\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, and $\mathrm{\Delta}y=5.3\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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) $\theta \cong \frac{5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}}{21\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}}\cong 12\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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% ($\text{speed}=1202\pm 5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{s}$).

## 3.4.2.

#### 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\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$; channel inner radius $\rho =1225\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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\times 2$ array [Fig. 6(b), $\mathrm{\Delta}x=14.8\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, $\mathrm{\Delta}y=4.6\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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\pm 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\pm 17\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{s}$, open squares, Fig. 6(c)] to 0.7% [$v=920\pm 7\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{s}$, filled squares, Fig. 6(c)].

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\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}\cong 0.7\text{\hspace{0.17em}}R$ from the channel axis), we find $v=408\pm 5\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{s}$ [Fig. 6(d), solid lines], that agrees well with the estimate $\cong {v}_{\mathrm{max}}[1-{(\frac{0.7\text{\hspace{0.17em}}R}{R})}^{2}]\cong 420\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{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\pm 3\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{s}$, not dissimilar from the speed in the field of view (b).

## 4.

## Discussion

## 4.1.

### 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 $\cong 50\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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\pm 10$ to $3\pm 0.7$ when the spot–spot distance decreases from 67 to $8\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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\times 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, ${\gamma}_{k,h\ne k}$ of the $\mathbf{\gamma}$ matrix and a geometric average of $k$’th and the $h$’th modes term $\sqrt{{\gamma}_{k,k}{\gamma}_{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\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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, $\zeta =\frac{1}{7}\sqrt{\sum _{n=5}^{11}{z}_{n}^{2}}$, does not increase dramatically when moving from the center of the FOV ($\zeta \cong 0.6$) to the periphery $\zeta \cong 0.7\pm 0.2$. Moreover, the most effective modes as measured by the mode amplitude averaged over the 9 ROIs, ${\sqrt{\u27e8{z}_{n}^{2}\u27e9}}_{\mathrm{ROI}}$, are $n=8$, 5, and 10 (${\sqrt{\u27e8{z}_{n}^{2}\u27e9}}_{\mathrm{ROI}}=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\pm 0.4\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ to $1.3\pm 0.2\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, Figs. 3(c)–3(e)] when we decrease the spot–spot distance in the grid from 67 to $8\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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 size^{56} 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.

## 4.2.

### 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\times 100\text{\hspace{0.17em}\hspace{0.17em}}\text{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\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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 $\text{voxel}\cong 350\times 350\times 1000\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{nm}}^{3}$ in a sample with a protein concentration $\cong 1\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{M}$ and a count per molecule $\cong 3000\text{\hspace{0.17em}\hspace{0.17em}}\text{counts}/\mathrm{s}$ (for a 10% collection efficiency), we should get $600\text{\hspace{0.17em}\hspace{0.17em}}\text{counts}/\mathrm{ms}$ or 240 electrons per pixel of an image acquired at 1 kHz (conversion factor $\cong 0.49\text{\hspace{0.17em}\hspace{0.17em}}\text{electrons}/\text{count}$ for Hamamatsu Orcad camera), corresponding to an acceptable level of noise $\cong 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\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{ms}$. For a $70\times 70\text{\hspace{0.17em}\hspace{0.17em}}\mu {\mathrm{m}}^{2}$ field of view, $0.14\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}\text{\hspace{0.17em}}\text{pixel}$ size, $512\times 512\text{\hspace{0.17em}\hspace{0.17em}}\text{pixels}$ (corresponding to about ${10}^{5}$ frames) the minimum acquisition time is then 9 min.

## 4.3.

### 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 vivo*^{8} 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\text{-}\mu \mathrm{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 $\sim 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 $\theta \cong \pm \frac{1\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}}{4.6\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}}\cong \pm 12\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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 $\cong \frac{4.6\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}}{14.8\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}}\cong 17\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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 $|\overrightarrow{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 number^{60} (i.e., large diffusion coefficients $D$ such that $|\overrightarrow{v}||\overrightarrow{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 $\overrightarrow{R}$, and therefore, $\alpha \ne 0$ [Eq. (14)]. A simple analysis of Eq. (14) suggests that the amplitude decay should follow an exponential law $\sim \mathrm{exp}(-{\alpha}^{2}/{\beta}^{2})$, where $\beta ={\omega}_{0}/|\overrightarrow{R}|$ is the angle under which the upstream spot sees the downstream spot (whose size is ${\omega}_{0}$). An analysis of the trend of the $\mathrm{ln}[{G}_{X12}({\tau}_{\mathrm{max}})]$ as a function of $\frac{1}{{\beta}^{2}}={|\overrightarrow{R}|}^{2}/{\omega}_{0}^{2}$ provides misalignment angles at most of $\alpha \cong 5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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 $\alpha $ value was found, indicating that the Peclet number in our case is large (in fact we estimate ${P}_{c}=|\overrightarrow{v}||\overrightarrow{R}|/D\cong \mathrm{19,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\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ and inner radius $\rho =1225\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$), we could have deviations from the axial flow such as Dean flow^{60} that could be detected in the image plane. However, since the Reynold number in our experiments is relatively low, ${R}_{y}\cong 0.3$, any secondary Dean flow^{60} in the radial direction is within the current uncertainty on the flow velocity direction. In fact, we expect^{45} a radial component of the velocity ${v}_{r}\cong {R}_{y}\frac{R}{\rho}|\overrightarrow{v}|$, or an angle of misalignment $\delta \theta \cong {R}_{y}\frac{R}{\rho}$ with respect to the laminar flow. In our case, this value $\delta \theta \cong 4.9\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{deg}$, well within the uncertainty due to the grid spacing, $\theta \cong \pm 12\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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 $\cong 1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{s}$. The spread in the CCF amplitude raises to 50% for a total acquisition time $\cong 0.4\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{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 $\delta {\tau}^{(0)}=2.5\pm 0.2$ to $\delta \tau =1.6\pm 0.1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{ms}$). The corresponding lagtime does not change appreciably with the application of the corrections but its uncertainty decreases, passing from ${\tau}_{\mathrm{max}}^{(0)}=9.7\pm 0.09\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{ms}$ to ${\tau}_{\mathrm{max}}=9.6\pm 0.03\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{ms}$. From these values, we can estimate the flow speed (the interspot distance is $d=14.7\pm 0.3\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$) that results ${v}^{(0)}=\frac{d}{{\tau}_{\mathrm{max}}^{(0)}}=1526\pm 14\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{s}$, when no aberrations correction is applied, and $v=1542\pm 6\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{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 ${\omega}_{0}$ through ${\omega}_{0}\cong v\delta \tau $ [see Eq. (14)]. The values of the spot size for the two cases are ${\omega}_{0}\cong 4.0\pm 0.3\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ and $2.5\pm 0.1\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{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\pm 0.2$ times in the range $120\le v\le 1200\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{s}$. Moreover, the accuracy increases threefold for the absolute speed value: for a nominal flow speed $v=120\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{s}$, we find $v=121.4\pm 0.06\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{s}$ (accuracy $\cong 1\%$) when the aberrations are corrected, and ${v}^{(0)}=124.0\pm 0.1\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}/\mathrm{s}$ (accuracy $\cong 3\%$) when no aberrations correction is applied [Fig. 5(e)].

## 5.

## Conclusions

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\text{\hspace{0.17em}\hspace{0.17em}}\frac{\mu \mathrm{m}}{\mathrm{s}}$, 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.

## Disclosures

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

## Acknowledgments

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.