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\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\mathrm{mm}$ and $f_2=140\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{pixel}$ resolution ($20\text{-}\mu \mathrm{m}\text{pixel}$ size) and 8-bit pixel depth. The SLM plane is conjugated through a $0.65\times$ beam reducer ($f_3=200\mathrm{mm}$ and $f_4=125\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,³⁵ while slightly displacing the position of the second lens of the beam reducer from the telescope condition.³¹

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 ($f_{\mathrm{LT}}=150\mathrm{mm}$) on a Hamamatsu CMOS Orca Flash 4.0 Camera (pixel size $6.5\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{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.³⁵ 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.³¹ The FL has the following expression:

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)³⁸ 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.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.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.³¹ 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:

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:

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

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 }_{\max }={\tau }_v$, for small values of the times ratio $\mathrm{\Gamma }={\tau }_v/{\tau }_D$. Numerical simulations proved⁶ that for $\mathrm{\Gamma }<0.1$ the discrepancy between ${\tau }_v$ and ${\tau }_{\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 }_{\max }\cong {\tau }_v\cos (\alpha )$. Moreover, the amplitude $G_{X12}({\tau }_{\max };{\tau }_D,{\tau }_v)$ of the CCF decreases as a function of the noncollinearity angle $\alpha$ as $G_{X12}({\tau }_{\max };{\tau }_D,{\tau }_v)\approx \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,⁵³ 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.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\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 50\mathrm{nm}\cong \frac{1}{7}\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.

Fig. 2

(a) Results of the spot identification algorithm on a two-photon image of a $5\times 5$ grid projected onto a Rhodamine slab (interspot spacing $21.1\times 15.8\mu \mathrm{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\times 5$ grid (interspot spacing $28.8\times 21.8\mu \mathrm{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).

3.2.

Aberrations Correction

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.

Fig. 3

(a) Juxtaposition of nine $5\times 5$ grids (aberrations corrected) projected onto a Rhodamine slab (two-photon fluorescence, interspot spacing $21.1\times 15.8\mu \mathrm{m}$ in the sample plane) with different combinations of $x$ shift ($-131.0$, 0, and $+131.0\mu \mathrm{m}$ in the sample plane) and $y$ shift ($-99.8,0$, and $+99.8\mu \mathrm{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\times 5$ grids projected onto a Rhodamine slab with increasing interspot spacing: (c) $7.7\times 6.1\mu \mathrm{m}$, (d) $34.8\times 26.3\mu \mathrm{m}$, and (e) $67.5\times 52.7\mu \mathrm{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).

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\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\mu \mathrm{m}$ to $2.67\pm 0.02\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\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\mu \mathrm{m}$ [Fig. 4(a)].

Fig. 4

(a) Single-axis profiles of spots obtained with aberration corrected grids of increasing number of spots $3\times 3$, $5\times 5$, and $7\times 7$, respectively, from bottom to top, each with interspots spacing equal to $27.8\times 20.7\mu \mathrm{m}$. The full-width at half-maximum of the spots measured in these three cases is $2.90\pm 0.03$, $2.75\pm 0.03$, and $2.33\pm 0.03\mu \mathrm{m}$, respectively. (b) The radial distribution for a noncorrected $5\times 5$ grid (same interspot spacing) has a full-width at $\text{half-maximum}=4.5\pm 0.8\mu \mathrm{m}$, almost double the nonaberrated case. (c) and (d) Two-photon images of a Convallaria sample obtained by employing a $13\times 13$ grid (spot–spot spacing $27.4\times 27.4\mu \mathrm{m}$) before aberrations correction [BC, (c)] and after the aberrations correction [AC, (d)] $\text{field of view}=356\times 356\mu \mathrm{m}$, excitation power $P=150\mathrm{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\times 5$ grid has been employed with a spacing $\mathrm{\Delta }x=14.0\mu \mathrm{m}$ and $\mathrm{\Delta }y=14.0\mu \mathrm{m}$ in the sample plane, shift steps $\mathrm{d}x=702\mathrm{nm}$ and $\mathrm{d}y=702\mathrm{nm}$ in the sample plane, corresponding to 2 CCD pixels in the CCD plane for a FOVs of $70.3\times 70.3\mu \mathrm{m}$: the excitation power was $P=150\mathrm{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).

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\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\mu \mathrm{m}$ in the sample plane) shifted with steps $\delta x=\delta y=702\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\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\mu \mathrm{m}$. High-resolution imaging of the same sample features with a Leica SP5 confocal microscope (resolution $330\pm 30\mathrm{nm}$) gives a width for these structures of $1.4\pm 0.5\mu \mathrm{m}$. These values indicate an optical resolution of the present setup of $1.5\pm 0.5\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\mu \mathrm{m}\times 800\mu \mathrm{m}$ size) and a horseshoe-shaped microchannel ($500\mu \mathrm{m}\times 500\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{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\mu \mathrm{m}$ and $\mathrm{\Delta }y=4.6\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)].

Fig. 5

CCF analysis on a portion of a straight microchannel. (a) Image (maximum signal projection of a time stack) of a $5\times 2$ grid (interspot spacing $14.8\times 4.6\mu \mathrm{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\pm 0.2$ and $1.6\pm 0.1\mathrm{ms}$ for the noncorrected and aberrations corrected cases, respectively. (c) Increase of the cross-correlation amplitude at the time-of-flight (${\tau }_{\max }$) as a function of the cut-off threshold (ROI size for the signal average, $3\times 3\text{pixels}$). The inset reports the CCFs for increasing values of the threshold. (d) Effect of the size ($m\times m$) of the ROI used to average the signal on the spot (threshold for cut-off = 25): $1\times 1$ (squares), $3\times 3$ (circles), $5\times 5$ (up-triangles), and $7\times 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 $G_{X12}({\tau }_{\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 ($\text{set values}=120$, 400, and $1200\mu \mathrm{m}/\mathrm{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\times 4$ grid (interspots spacing $7.0\times 5.3\mu \mathrm{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\pm 5\mu \mathrm{m}/\mathrm{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 $|\overrightarrow{R}|$ normalized to the size of the spot ${\omega }_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 ($\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{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{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{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\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\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\mu \mathrm{m}$, and $\mathrm{\Delta }y=5.3\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\mu \mathrm{m}}{21\mu \mathrm{m}}\cong 12\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\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\mu \mathrm{m}$; channel inner radius $\rho =1225\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\mu \mathrm{m}$, $\mathrm{\Delta }y=4.6\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\mu \mathrm{m}/\mathrm{s}$, open squares, Fig. 6(c)] to 0.7% [$v=920\pm 7\mu \mathrm{m}/\mathrm{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\times 2$ grid (interspot spacing $14.8\times 4.6\mu \mathrm{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 ($\text{frame rate}=855\mathrm{fps}$). The flow speed recovered from Gaussian fit was $920\pm 7\mu \mathrm{m}/\mathrm{s}$ and $900\pm 17\mu \mathrm{m}/\mathrm{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\pm 5\mu \mathrm{m}/\mathrm{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) ($\text{frame rate}=986\mathrm{fps}$). The speed retrieved by the best Gaussian fit (solid line) of the pair H–C of spots was $457\pm 3\mu \mathrm{m}/\mathrm{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\mu \mathrm{m}\cong 0.7R$ from the channel axis), we find $v=408\pm 5\mu \mathrm{m}/\mathrm{s}$ [Fig. 6(d), solid lines], that agrees well with the estimate $\cong v_{\max }[1-{(\frac{0.7R}{R})}^2]\cong 420\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\mu \mathrm{m}/\mathrm{s}$, not dissimilar from the speed in the field of view (b).

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\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\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\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{⟨z_n^2⟩}}_{\mathrm{ROI}}$, are $n=8$, 5, and 10 (${\sqrt{⟨z_n^2⟩}}_{\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\mu \mathrm{m}$ to $1.3\pm 0.2\mu \mathrm{m}$, Figs. 3(c)–3(e)] when we decrease the spot–spot distance in the grid from 67 to $8\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⁵⁶ 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.³¹ The total acquisition time for the image shown in Fig. 4(e) ($100\times 100\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\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{\mathrm{nm}}^3$ in a sample with a protein concentration $\cong 1\mu \mathrm{M}$ and a count per molecule $\cong 3000\text{counts}/\mathrm{s}$ (for a 10% collection efficiency), we should get $600\text{counts}/\mathrm{ms}$ or 240 electrons per pixel of an image acquired at 1 kHz (conversion factor $\cong 0.49\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\mathrm{ms}$. For a $70\times 70\mu {\mathrm{m}}^2$ field of view, $0.14\mu \mathrm{m}\text{pixel}$ size, $512\times 512\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,⁵⁹ or in vivo⁸ 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.³¹ 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\mu \mathrm{m}}{4.6\mu \mathrm{m}}\cong \pm 12\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\mu \mathrm{m}}{14.8\mu \mathrm{m}}\cong 17\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⁶⁰ (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 \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 $\ln [G_{X12}({\tau }_{\max })]$ as a function of $\frac{1}{{\beta }^2}={|\overrightarrow{R}|}^2/{\omega }_0^2$ provides misalignment angles at most of $\alpha \cong 5\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\mu \mathrm{m}$ and inner radius $\rho =1225\mu \mathrm{m}$), we could have deviations from the axial flow such as Dean flow⁶⁰ 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⁶⁰ in the radial direction is within the current uncertainty on the flow velocity direction. In fact, we expect⁴⁵ 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\mathrm{deg}$, well within the uncertainty due to the grid spacing, $\theta \cong \pm 12\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\mathrm{s}$. The spread in the CCF amplitude raises to 50% for a total acquisition time $\cong 0.4\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\mathrm{ms}$). The corresponding lagtime does not change appreciably with the application of the corrections but its uncertainty decreases, passing from ${\tau }_{\max }^{(0)}=9.7\pm 0.09\mathrm{ms}$ to ${\tau }_{\max }=9.6\pm 0.03\mathrm{ms}$. From these values, we can estimate the flow speed (the interspot distance is $d=14.7\pm 0.3\mu \mathrm{m}$) that results $v^{(0)}=\frac{d}{{\tau }_{\max }^{(0)}}=1526\pm 14\mu \mathrm{m}/\mathrm{s}$, when no aberrations correction is applied, and $v=1542\pm 6\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\mu \mathrm{m}$ and $2.5\pm 0.1\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\mu \mathrm{m}/\mathrm{s}$. Moreover, the accuracy increases threefold for the absolute speed value: for a nominal flow speed $v=120\mu \mathrm{m}/\mathrm{s}$, we find $v=121.4\pm 0.06\mu \mathrm{m}/\mathrm{s}$ (accuracy $\cong 1\%$) when the aberrations are corrected, and $v^{(0)}=124.0\pm 0.1\mu \mathrm{m}/\mathrm{s}$ (accuracy $\cong 3\%$) when no aberrations correction is applied [Fig. 5(e)].