Quantitative Lateral and Axial Flow Imaging with Optical Coherence Microscopy and Tomography References and Links

Optical coherence tomography (OCT) and optical coherence microscopy (OCM) allow the acquisition of quantitative three-dimensional axial flow by estimating the Doppler shift caused by moving scatter-ers. Measuring the velocity of red blood cells is currently the principal application of these methods. In many biological tissues, blood flow is often perpendicular to the optical axis, creating the need for a quantitative measurement of lateral flow. Previous work has shown that lateral flow can be measured from the Doppler bandwidth, albeit only for simplified optical systems. In this work, we present a generalized model to analyze the influence of relevant OCT/OCM system parameters such as light source spectrum, numerical aperture and beam geometry on the Doppler spectrum. Our analysis results in a general framework relating the mean and variance of the Doppler frequency to the axial and lateral flow velocity components. Based on this model, we present an optimized acquisition protocol and algorithm to reconstruct quantitative measurements of lateral and axial flow from the Doppler spectrum for any given OCT/OCM system. To validate this approach, Doppler spectrum analysis is employed to quantitatively measure flow in a capillary with both extended focus OCM and OCT. Optical coherence tomography-principles and applications, " Rep. Mobility and transverse flow visualization using phase variance contrast with spectral domain optical coherence tomography, " Opt. Speckle variance detection of microvasculature using swept-source optical coherence tomography, " Opt. In vivo volumetric imaging of vascular perfusion within human retina and choroids with optical micro-angiography, " Opt.ric angiography of cortical microvasculature with optical coherence tomography, " Opt. Extended focus high-speed swept source OCT with self-reconstructive illumination, " Opt. In vivo bidirectional color Doppler flow imaging of picoliter blood volumes using optical coherence tomography, " Opt. Optical Doppler tomographic imaging of fluid flow velocity in highly scattering media, " Opt. " Real-time measurement of in vitro flow by Fourier-domain color Doppler optical coherence tomography, " Opt. Resonant Doppler flow imaging and optical vivisection of retinal blood vessels, " Opt. Flow velocity estimation using joint spectral and time-domain optical coherence tomography, " Opt. Doppler optical micro-angiography for volumetric imaging of vascular perfusion in vivo, " Opt. Quantitative cerebral blood flow with optical coherence tomography, " Opt. retinal blood flow measurement with ultrahigh speed swept source/Fourier domain OCT, " Biomed. of microvasculature by dual-beam phase-resolved Doppler optical coherence tomography, " Opt. Three-dimensional quantitative imaging of retinal and choroidal blood flow velocity using joint spectral and …


Introduction
Optical coherence tomography (OCT) is an interferometric imaging technique that allows threedimensional label-free imaging of biological tissues and cells [1].Besides structural imaging, OCT has been used extensively for functional imaging, where blood flow imaging in living tissue is arguably one of the most important applications.Blood flow imaging methods based on OCT can be divided into two categories: angiographic methods that image the structure of the vasculature in tissue [2][3][4][5][6][7][8], and velocimetric methods measuring the velocity of red blood cells [9][10][11][12][13][14][15][16][17][18].While the former techniques use some form of motion contrast to yield a qualitative image of blood vessels, the latter allow a quantitative estimation of the blood flow speed.Most quantitative methods are ultimately based on measuring the Doppler frequency shift of light scattered by moving particles, but use different approaches to estimate the speed of the scatterers from the acquired data.All these methods require the acquisition of a time-series of A-scans or B-scans.Hence, the highest Doppler frequency that can be measured without aliasing is given by half the acquisition rate.The most straightforward approach is to evaluate the local phase difference between time samples, from which the Doppler frequency is then obtained.In joint spectral and time domain OCT, Fourier transformations are performed along both time and wavenumber dimensions, yielding the spatially resolved Doppler spectrum [14,19,20].At each location, the Doppler shift is estimated as the frequency corresponding to the maximum magnitude position or to the center of mass using circular statistics [21].Autocorrelation-based methods involve a similar acquisition scheme but calculate the Doppler shift by searching the maximum in the power spectral density calculated as the Fourier transform of the autocorrelation [16].
In all these quantitative methods, the estimated Doppler frequency is eventually used to calculate the axial velocity component of moving scatterers.However, in a number of biological applications (e.g.neurology and ophthalmology), many capillaries are oriented nearly perpendicular to the optical axis, making the axial flow component too small to be measured in small vessels such as capillaries.It is therefore crucial to provide a reliable measurement of the lateral component.
The lateral velocity component can be estimated from the width of the Doppler frequency spectrum [22][23][24][24][25][26][27].Indeed, because wave vectors in the beam illuminating the sample and in the collected backscattered light have an angular spread with respect to the optical axis, they will experience different Doppler shifts, causing broadening of the Doppler spectrum.This effect is called spectral broadening due to beam geometry [28].Proskurin et al. argued that the finite transit time of scatterers through the focal volume causes spectral broadening [25].This effect can be shown to be equivalent to broadening due to beam geometry [28].The equivalence can be understood intuitively by considering the effects of increasing the numerical aperture (NA) of an optical system.A higher NA results in a larger angular spread of wave vectors, increasing broadening of the Doppler spectrum due to the beam geometry argument.Equivalently, the increased angular spread of wave vectors yields an improved lateral resolution, thereby reducing transit time through the focus, which results in broadening of the Doppler spectrum.
A conceptually different method to measure axial and lateral flow velocity components as well as diffusion was presented by Lee et al. [44].Their method is based on dynamic light scattering measurements with OCT and involves fitting the autocorrelation of the OCT signal to a model.While this method allows to extract a large amount of information about the local flow properties, the non-linear fitting of the model's many parameters requires an intricate algorithm resulting in long computation times.
Phase noise in the interferometer creates Doppler broadening even in the absence of flow, thus imposing a limit on the minimum detectable lateral flow velocity [14,29].Therefore, the larger the Doppler frequency broadening for a given flow velocity vector, the better the sensitivity to small flow changes.Higher NA increases Doppler broadening due to the lateral flow component [23].Hence, to image slow lateral flow or to detect small flow changes, a higher NA is advantageous.In this case, optical coherence microscopy (OCM) can have an advantage over OCT.For OCM systems employing Gaussian beams however, this results in a reduced depth of focus.In this work, we use extended focus optical coherence microscopy (xfOCM) to counter this problem.xfOCM uses a Bessel mode to illuminate the sample and a Gaussian mode to collect scattered light [30].The Bessel-Gauss configuration achieves a high lateral resolution of 1.3 µm, maintained over an extended depth of focus of 400 µm.
To date, the models developed to analyze Doppler broadening in OCT have all assumed Gaussian modes and light source spectra.In this work, we develop a model that allows arbitrary light source spectra, and arbitrary illumination and detection modes.This allows modeling of arbitrary axial and lateral parts of the point spread function, respectively.Given the optical system's parameters, our model can be used to correctly combine the measured central Doppler frequency and Doppler broadening to provide quantitative images of axial and lateral flow velocities.To validate the model, we perform quantitative imaging of lateral and axial flow in a capillary channel using both an OCT system with Gaussian modes and an xfOCM instrument.

Fourier domain coherent imaging for static scattering
Typical Fourier domain implementations of OCT record an interferogram as a function of wavenumber at each position in a lateral raster scan.After background subtraction and assuming the backscattered light intensity is small compared to the reference light intensity, the signal of interest in each wavenumber channel is where U r is the reference field amplitude and U m the scattered field amplitude measured by the system.Each wavenumber channel can be treated as being monochromatic.The optical system  We will assume that the objective lens obeys Abbe's sine condition [31].For scalar waves and under the Debye approximation [32], the incident light field in the focal volume at a position r measured from the focal point, can be written as a superposition of plane waves with different angular directions [33]: where f is the focal length of the objective, k i (p i ) describes the illumination wave vectors at the exit principal plane and S(k) is the spectral envelope of the field.In the notation U in (r; k), r is a variable for the function U in whereas k is a parameter indicating the (monochromatic) wavenumber channel.The direction cosines of k i are functions of p i = (a i , b i ): The plane waves are scattered upon interaction with the sample.Using the first order Born approximation [34], the scattered field in the principal plane is where r d (p d ) describes the exit principal plane of the objective (as shown in Fig. 1(a)) and F(r; k) is the scattering potential.The scattering potential can be expressed in terms of the sample susceptibility (assumed to be non-dispersive): Because Finally, the coupling into the detection mode needs to be considered: For a detailed analysis of the tomogram reconstruction from Eq. ( 7), we refer to the works of Villiger et al. [35] and Sheppard et al. [36].In the case of NA → 0, both k d and k i are parallel to the optical axis, as shown in Fig. 1(b).Then, k d = −k i = k1 z and hence the measured scattered field would be: If the scattering potential is taken to be a set of scattering layers, the classical layer model for OCT image formation is retrieved.In this work however, we will start from Eq. ( 7) without further approximation because the NA is crucial to measure lateral and axial flow components.

Fourier domain coherent imaging for dynamic scattering
When scatterers in the sample are moving, the scattering potential becomes a function of time and can be separated in a static and dynamic part: The dependence of the dynamic term F d on the time parameter t is assumed to be on a much slower time-scale than the period of the light wave, such that the monochromatic approximation leading to Eq. (4) remains valid.A sample containing different ensembles of uniformly moving rigid scatterers can be modeled by: where the v n are the velocity vectors of the ensembles and n runs over the contributing ensembles.To evaluate the volume integral in Eq. ( 6) for each of the ensembles of moving scatterers, we can adopt a coordinate system r 0 = r − v n t that is moving along with the ensemble.The time-dependent part U m (t; k) of the measured field is then: In Eq. (11), the Fourier transformed scattering potential can be identified: From the measurement of the field U m , OCT systems attempt to reconstruct the depth resolved scattering potential of the sample at each position in a lateral raster scan.Typically, an inverse Fourier transform of Eq. ( 1) is calculated over k to this end.For this reconstruction to be accurate at some distance from the focal point, either the range of K(p d , p i ) must be limited by m d (p d ) and m i (p i ) such that r 0 • K effectively measures the optical path length along the depth direction (i.e.r 0 • K ≈ 2zk), or advanced reconstruction algorithms need to be applied.The former solution is valid within good approximation in the focal volume of low NA or extended focus systems [35].The latter solution has been investigated under the name of interferometric synthetic aperture microscopy (ISAM) [37][38][39].In this work, we will assume either one of the solutions is adopted to ensure correct reconstruction.Defocussing effects will therefore not be considered.
Finally, the time-dependent part of the measured field can be rewritten as: In this expression, the functions m d (p d ) and m i (p i ) can be interpreted as distribution functions for the variables p d and p i .Integration over these variables creates the Doppler spectrum as follows.For each realization of p d and p i , the corresponding illumination and detection wave vectors can be found, as shown in Fig. 1(a).Their direction cosines are given by Eq. ( 3).Every such pair of illumination and detection wave vectors determines a Doppler frequency shift per wavenumber channel: The weight of each frequency in the Doppler spectrum is then given by: In the extreme case of zero NA shown in Fig. 1(b), we find that the distribution reduces to a single Doppler frequency, as expected.Moreover, lateral flow will not be detected in that case, since v n • K = 0. Conversely, when the NA is increased, K spans a range of angles due to the spread of wave vectors in the illumination and detection modes, as shown in Fig. 1(c).Therefore, a broadened Doppler frequency distribution will be seen and lateral flow can be measured.An analysis of the Doppler spectrum allows the measurement of both flow components, as will now be shown.

Doppler spectrum per wavenumber channel: mean and standard deviation
To analyze the Doppler spectrum per wavenumber channel, the weighted mean µ[ f D ; k] and standard deviation σ [ f D ; k] of the Doppler frequency are evaluated: For the sake of simplicity, a single point scatterer has been taken into account ( Fd,n In other words, back-scattering is approximated to be uniform within the system's aperture.The model in principle also allows different types of scatterers by incorporating Fd,n (K) in the Doppler frequency weight, but such an analysis is beyond the scope of this work.Since most optical systems exhibit rotationally symmetric illumination and detection modes, polar coordinates (ρ, θ ) will be used to calculate the integrals in the entrance principal plane.Considering one ensemble of point scatterers with velocity vector v = (v x , v y , v z ), we find for the mean Doppler frequency for one wavenumber channel: where n is the refractive index of the medium and the following notation for integrals of the illumination mode has been adopted: with R the radius of the entrance pupil.Analogous notation is used for integrals of the detection mode.Eq. ( 18) confirms that the mean Doppler frequency is determined solely by the axial flow component.Moreover, when a moderate NA objective is used, m(ρ) is significantly different from zero only for ρ f and therefore M f ≈ f M 1 .This approximation is further verified in Fig. 2(b).The mean Doppler frequency then reduces to the known formula: In general, the variance of the Doppler frequency is dependent on all flow components: with The magnitudes of the ratios M f /M 1 and M 3 /M 1 , are shown in Fig. 2 for Gaussian and Bessel modes.The dependence of the Doppler broadening on A, f and B can be interpreted as follows.
As mentioned in the introduction, it has been shown that broadening of the Doppler spectrum due to beam geometry is equivalent to broadening due to the finite transit time of scatterers moving through the focal volume [28].Hence, Doppler broadening due to lateral flow is higher when the lateral resolution is higher.Higher lateral resolution requires increasing the NA, which in turn will increase A, as shown in Fig. 2(c).This explains the presence of A in the first term of Eq. ( 23).The second term measures broadening due to pure axial flow, which is caused by the angular spread of wave vectors with respect to the axial component of the velocity vector.Wave vectors with different angles measure different Doppler frequencies.Equivalently, this can be interpreted by the finite depth of focus of the system, which limits the focal volume axially.The depth of focus therefore creates Doppler broadening due to the finite axial transit time of scatterers moving through the focal volume.For high NA Gaussian systems, the depth of focus will be reduced and the factor C increases, as shown in Fig. 2(d).However, this effect is only of importance when the flow is mostly axial, since C A/2, as can be seen by comparing Fig. 2(c) with (d).Note that for extended focus systems with ideal Bessel illumination and detection modes (m i,d (ρ) = δ (ρ − ρ 0 )), we find C = 0 for all NA.
Figure 3 shows the mean and standard deviation of the Doppler frequency calculated using Eq. ( 18) and Eq. ( 23) for different angles of the velocity vector, and for different NAs.The mean Doppler frequency shown in Fig. 3(c) is seen to be almost independent of the NA, as predicted by the approximation leading to Eq. (22).For pure axial flow, however, we see that the mean Doppler frequency is slightly lower for higher NA.This is due to the smaller mean z-component of K for higher NA.On the other hand, a higher NA increases the Doppler broadening, as evidenced in Fig. 3(d).When a typical OCT system of NA = 0.05 is compared to a system  with NA = 0.3, the Doppler broadening for lateral flow is increased about sixfold.Fig. 3(d) confirms that even for pure axial flow, some Doppler frequency broadening remains, especially for higher NA.This is caused by the second term of Eq. ( 23), as explained above.

Spatially resolved Doppler spectrum
To yield depth-resolved imaging of the Doppler spectrum, we build on the method of joint spectral and time domain OCT (jstdOCT) [14].After subtraction of the background spectrum, the interferograms recorded by the spectrometer are first mapped from the wavelength to the wavenumber space.Then, the inverse Fourier transform is calculated along the wavenumber axis: where 2l represents optical path length, and the factor 2 takes into account the double path due to the reflection configuration.As discussed above, we will assume that Ũm (t, l) correctly represents the depth resolved sample structure.The Doppler spectrum can be estimated from the Fourier transform of Eq. (25) along t.After this Fourier transform, we have a measurement of the Doppler spectrum for each depth resolution.Indeed, the integral over k in Eq. ( 25) selects one resolution volume at depth l.More details on the full algorithm are given in section 3.4.
The integration over the k-spectrum will cause the Doppler frequency distribution to broaden.Indeed, f D in Eq. ( 14) depends linearly on k, the scattering potential in the Doppler frequency weight in Eq. ( 15) contains the factor k 2 and upon inspection of Eq. ( 13), it is seen that the additional weighting factor kS(k) must be taken into account.The effect on the weighted mean and standard deviation of the Doppler frequency for a point scatterer as calculated in section 2.3, can now be evaluated.
The mean Doppler frequency can be calculated using the law of total expectation: In the case of point scatterers, we can substitute Eq. ( 18) for µ[ As before, for moderate NA, M f ≈ f M 1 .When S(k) is a Gaussian function with mean k 0 and standard deviation k σ , the moments of S(k) can easily be calculated.Moreover, k 2 σ k 2 0 for typical OCT/OCM spectra, which allows to simplify the ratio of the fourth to the third moment of S(k) to k 0 .With these two approximations, Eq. ( 27) reduces again to nv z k 0 /π.
Similarly, the variance is found using the law of total variance: Substituting Eq. ( 23) for To interpret the different terms of Eq. ( 29), we first assume S(k) to be Gaussian with central wavenumber k 0 , standard deviation k σ and k 2 σ k 2 0 , which is valid for typical OCT/OCM spectra.Finally, with M f ≈ f M 1 for moderate NA we find: The two terms of Eq. ( 30) can be interpreted as follows.The first term T 1 is the standard deviation for a single wavenumber channel at the central wavenumber as given by Eq. ( 23); it accounts for broadening due to the angular spread of incident and scattered wave vectors, as discussed above.Note that T 1 is independent of spectral width.The second term T 2 measures an extra broadening due to the axial flow component, and depends on the spectral width k σ .It can be interpreted using the transit time argument: a larger spectral bandwidth creates a higher axial resolution, resulting in a shorter axial transit time, and therefore a larger Doppler broadening.Indeed, the axial resolution is given by the coherence length l c = 1/(nk σ ).
Figure 4 shows the mean and standard deviation (from Eq. ( 29)) of the Doppler frequency for different NA and light source bandwidth.The mean Doppler frequency is independent of these two parameters, as seen in Fig. 4(a) and (c), and as predicted by the approximations.The behavior of the standard deviation in Fig. 4(b) and (d) can be understood in terms of spatial resolution and the transit time argument.When axial resolution (determined by the source spectrum) and lateral resolution (determined by the NA) are equal, the Doppler frequency broadening is equal for axial and lateral flow components, since the transit time through the focal volume is then independent of flow direction.Conversely, an asymmetric focal volume will have higher Doppler broadening for either axial flow (better axial than lateral resolution) or for lateral flow (better lateral than axial resolution).

Approximations for Gaussian modes and low NA
To compare our formula in Eq. ( 30) to previous work [23,24,27], we need to set the illumination and detection modes to narrow (low NA) and equal Gaussians with waist w 0 R at the principal plane.Under this condition, A = 2w 2 0 and C ≈ 0. The beam waist after focussing by the objective is then w 0 = 2 f /(nk 0 w 0 ), yielding for the Doppler broadening: where 1/2 is the transverse flow component.The first term in Eq. ( 31) corresponds to the formula for Doppler broadening as given by Ren et al. [23] and Piao etal.[24].Their works do not include broadening due to the axial flow component.Srinivasan et al. calculate the broadening from the field autocorrelation for equal axial and lateral resolution [27]; their result is equivalent to Eq. (31).In conclusion, with the right approximations for low NA OCT, our model reduces to the results in previously published work.Table 1 summarizes the different approximations for the mean and standard deviation of the Doppler spectrum.

Doppler spectrum: numerical simulation
To investigate the shape of the Doppler spectrum, a numerical simulation is performed.In the entrance principal plane, two sets of uniformly distributed random coordinates are chosen: p i for the illumination mode and p d for the detection mode.For each combination of p i and p d , the Doppler frequency and its weight are calculated using Eq. ( 14) and Eq. ( 15), including the spectral weighting factor for point scatterers k 3 S(k) for a Gaussian spectrum.This calculation is repeated for each wavenumber channel.The resulting data are then binned according to their frequency and represented in a weighted histogram.

xfOCM
The xfOCM set-up is based on a Mach Zehnder interferometer and has decoupled illumination and detection paths.Light from a broad-bandwidth source (Ti-Sa laser, Femtolasers, λ 0 = 780 nm, ∆λ = 120 nm) is collimated and split into reference and sample beams.In the sample arm, the beam passes through an axicon lens and a telescope.The resulting Bessel-like illumination beam is guided through a relayed beam scanning system.The scanned beam enters a microscope stage, consisting of a tube lens and a 10x Zeiss Neofluar objective with a NA of 0.3.The Bessel-like illumination beam provides a uniform, high lateral resolution (1.3 µm) over an extended depth of field (400 µm).The axial resolution in tissue is 2.5 µm.The system operates with a Bessel illumination mode and a Gaussian detection mode, as shown in Fig. 2(a).A custom-built spectrometer with a high-speed line detector (Basler Sprint spL4096-140km) records the interference spectrum.For the flow measurements presented here, an A-scan rate of 20 kHz was used.

OCT
The OCT system is based on a Michelson interferometer with equal Gaussian illumination and detection modes and employs a broad-bandwidth source (Ti-Sa laser, Femtolasers, λ 0 = 780 nm, ∆λ = 120 nm).The sample objective lens has a focal length of f = 25 mm.The Gaussian mode has a full width at half maximum of 2.4 mm in the back focal plane of this lens, creating a lateral resolution of 10 µm.A custom-built spectrometer with a high-speed line detector (Basler Sprint spL4096-140km) records the interference spectrum.For the flow measurements, an A-scan rate of 100 kHz was used.More details about this set-up can be found in a previous publication [40].

Flow system
The flow system consists of a glass capillary tube with a measured cross-section of 4.95 × 10 4 µm 2 for xfOCM and 3.14 × 10 4 µm 2 for OCT experiments.A syringe pump allows to set a constant flow speed through the capillary.The system is filled with a solution of polystyrene beads (diameter 0.2 µm) in water for the xfOCM-based measurements and with an intralipid solution for the OCT-based measurements.

Data acquisition & algorithm
A lateral raster scan is performed with oversampling (xfOCM: 16× Nyquist frequency; OCT: 125× Nyquist frequency) along the fast scan axis.Due to the oversampling, repeated A-scans are taken within the same lateral resolution volume such that they can be seen as separated in time, not space.The acquired interference spectra I (λ , x(t)) from each lateral resolution volume are first windowed along the time-dimension using a Hanning window, yielding I w (λ , x,t).Then, following the joint spectral and time-domain method [14], the spectral background is subtracted, the spectra are mapped to the wavenumber domain and the absolute value of the two-dimensional fast Fourier transform (2D-FFT) is calculated, yielding the Doppler frequency distribution as a function of depth.Figure 6 summarizes the data acquisition protocol and algorithm.For the xfOCM measurements a time-window size of 32 points is used, for OCT 256 points.To estimate the mean Doppler frequency, the center of mass of the Doppler spectrum is calculated using circular statistics [21].The Doppler frequency standard deviation can be estimated by calculating the weighted standard deviation directly, but we have found the following fitting procedure to be more reliable in the presence of noise.First, to prevent aliased Doppler frequencies from biassing the fit, each Doppler spectrum is circularly shifted to have its center of mass at 0 Hz.Then, the Doppler spectrum is fitted with a non-linear least squares algorithm (Matlab, The Mathworks, Inc.) to This empirical function was found to produce a good fit for both OCT and xfOCM Doppler frequency distributions, as shown in Fig. 6(e) for typical xfOCM data.The fitting parameter c is included to fit the (broader) xfOCM distributions; for OCT we use c → ∞.Indeed, the numerical simulations in Fig. 5(c) and Fig. 5(d) show that Doppler spectra in xfOCM are broader than in OCT.The Doppler frequency standard deviation is then calculated using p( f D ) as the weighting function.Lateral and axial flow components can then be found from the mean and standard deviation of the Doppler frequency using Eq. ( 27) and Eq. ( 29).To calculate the system-related parameters (i.e. the spectrum, illumination and detection modes) needed in these equations, we used the design parameters of the systems.A Matlab implementation of the above algorithm is provided, see the supporting information for details.Note that the same principles can be applied along the slow scan axis, by acquiring repeated B-scans; see Grulkowski et al. for a detailed description of scan protocols [20].Eq. ( 29) v t Step 2 I w(λ,t) Step 1 I(λ,x(t)) Step Step 5 v z (z) v t (z) Eq. ( 27) Eq. ( 29)  27) and Eq. ( 29), or their appropriate approximations.

Experimental results
To validate the model and algorithm, we performed flow measurements at different flow rates set by a syringe pump through a glass capillary using two optical set-ups: xfOCM and OCT.
The lateral scanning was perpendicular to the flow direction, such that its influence on Doppler broadening could be removed using v 2 t = v2 t − v 2 s , where v t is the true transverse flow velocity, vt the transverse velocity measured from the Doppler broadening and v s the scan speed.Another solution would be to make v s v t , or to stop the scan at each lateral position.As can be seen in 7(a), the capillary is somewhat asymmetric, with an elliptic cross-section.The axial diameter is 200 µm, but the lateral diameter is larger, accounting for a larger cross-section.The intensity gradient over the depth of the flow channel seen in 7(a) is due to three reasons.Firstly, although the extended focus maintains the lateral resolution over 400 µm, it does not have constant intensity over the whole depth, a well-known property of Bessel beams [42].Secondly, the Gaussian detection mode employed by xfOCM   creates and additional intensity gradient over depth.Thirdly, scattering of the illumination beam by the sample reduces the intensity with depth.

xfOCM
For each lateral scan position, µ[ f D ] and σ [ f D ] are converted to lateral and axial flow profiles as a function of depth, as shown in Fig. 8(a) and (b) for the center of the capillary.By assuming laminar flow, the volume flow rate F set by the syringe pump is converted to flow velocities in the capillary.The maximum velocity v max of the parabolic flow profile in a cylindrical capillary is then found from v max = 2F/S, with S the measured cross-section of the capillary.The flow velocities set by the syringe pump can now be compared to the velocities measured using our model and algorithm, as shown in Fig. 8(c) and (d) for an angle α = 81 • .The measurements of both axial and lateral flow components are seen to be consistent with the values expected from the flow system parameters.The measured flow speeds in Fig. 8(a) and (b) show some deviation from the parabolic flow profiles most likely due to experimental noise and/or deviations from the supposed laminar flow profile.Also, note that when a particular axial velocity is higher than the expected value, the corresponding transverse velocity component is higher as well, indicating that these errors are most probably due to an inaccuracy in the actual flow rates set by the syringe pump.

OCT
Figure 9(a) and (b) shows the depth profiles measured with the OCT system on a capillary with angle α = 87 • .Due to the lower lateral resolution of this system (i.e.smaller Doppler broadening), it is not as sensitive to slow lateral flow velocities as xfOCM.Hence, the flow rates used for these measurements were about ten times higher and the flow angle smaller to create larger lateral flow velocities.Nevertheless, the measurements of both axial and lateral flow components in Fig. 9(c) and (d) are again seen to match the values expected from the flow system parameters.

Discussion & Conclusion
We have introduced a general model for the Doppler frequency spectrum caused by moving scatterers that is applicable to any OCT or OCM system, given its parameters.Once all system related factors have been calculated (this must be done only once per set-up), Eq. ( 27) and Eq. ( 29) (or their appropriate approximations, see table 1) yield a straightforward relation between mean and standard deviation of the Doppler spectrum and the flow velocity components.
For the specific case of low NA OCT and Gaussian source spectra, our model yields the same equations as previously published models [23,24,27].The model presented here therefore provides a generalization of these earlier models to arbitrary point spread functions, through the possibility of arbitrary illumination and detection modes, and light source spectra, respectively.This generalization is important for high NA OCM systems as well as xfOCM and dark-field OCM [45].Additionally, the model does not make assumptions about the flow angle.
While point scatterers were investigated here, the effect of different types of scatterers could be investigated by including the scattering potential Fd,n (K) as a weight in the calculation of the mean and standard deviation of the Doppler spectrum.In this way, the Doppler spectrum created by moving erythrocytes in blood could be investigated.
The model assumes that the scattering structures are rigid bodies, moving with a constant velocity during the acquisition time.This assumption excludes effects such as Brownian motion or significant velocity gradients within a resolution volume.For the model to be valid, these effects must be small in comparison to the average velocity.
To measure the Doppler spectrum, we presented an algorithm building on the joint spectral and time domain method, thereby adopting the advantages of this method in terms of precision and signal-to-noise ratio over phase-resolved methods [14].To estimate the mean frequency of the Doppler spectrum, either circular statistics or the standard definitions can be used [21].A fitting procedure was used to allow a reliable estimation of the Doppler frequency standard deviation.For an interesting comparison of Doppler frequency estimators in the presence of noise, we refer to the recent work by Chan et al. [41].An interesting subject for future research would be to evaluate and compare the performance of OCT/OCM systems in terms of quantitative flow imaging.To this end, the influence of noise should be included in the model.
Our model and algorithm were successfully applied to image axial and lateral flow in a capillary using two different system configurations.No system modifications (digital or optical) to the point spread function are needed to use our method, only the scan protocol needs to be adapted for over-sampling of the lateral resolution.As mentioned before, the scan protocol can be modified to perform Doppler spectral analysis along the slow scan axis, allowing to concentrate on a lower set of velocities [20,27].By allowing arbitrary pupil functions, this work now enables high resolution extended-focus systems to be used for quantitative imaging of both axial and lateral flow components.Especially the ability to quantify the lateral flow component will be instrumental in a number of biological applications, such as retina or brain imaging where blood vessels are often oriented perpendicular to the optical axis.Furthermore, xfOCM has been shown to achieve a resolution and penetration depth comparable to multi- photon microscopy but with significantly faster acquisition rates [43].The additional capability to quantify blood flow velocities will allow the investigation of the three-dimensional vascular functionality at speeds currently unattainable by multiphoton microscopy.Our future work will therefore be focused on the quantitative imaging of lateral and axial flow in biological tissue using the concepts developed here.

Supporting information
A Matlab implementation of the algorithm developed here together with example data can be downloaded from our website at http://lob.epfl.ch/under "Research/Quantitative flow imaging".
is characterized by its (achromatic) illumination mode m i (a, b) and detection mode m d (a, b), parametrized in the entrance principal plane with coordinates p = (a, b), as in Fig. 1(a).

Fig. 1 .
Fig. 1.Schematic representations of wave vectors and the principal planes of an objective.(a) The axes and definitions used in the model.P and Q are the entrance and exit principal planes of the objective, respectively.(b) A low NA system has a small angular spread of wave vectors; lateral flow is therefore not detected because the flow velocity vector v n is then perpendicular to K = k i (p i ) − k d (p d ).Equal Gaussian illumination and detection modes are shown in the schematic.(c) Higher NA systems have a larger angular spread of wave vectors and the scalar product is no longer always zero, allowing the detection of lateral flow.

Fig. 2 .
Fig. 2. (a) Illumination or detection mode in the principal plane for Gaussian mode (blue) and Bessel mode (green) in a system with NA = 0.3.(b) The ratio M f /M 1 normalized to the focal length f for Gaussian and Bessel modes as a function of NA.For lower NA, M f /M 1 approach f .(c) M 3 /M 1 normalized to f 2 as a function of NA.(d) C = 2 f 2 − A − B normalized to f 2 as a function of NA for Gaussian illumination and detection modes.Note that C M3/M1.

Fig. 3 .
Fig. 3. (a) The velocity vector v of a point scatterer makes an angle α with the optical axis.(b) Equal Gaussian illumination and detection modes for an optical system with f = 16.4 mm.NA = 0.05, 0.3 and 0.5 for the blue, green and red curves respectively.(c) Mean Doppler frequency as a function of α.Parameters: λ = 780 nm, |v| = 1 mm/s, n = 1.33.(d) Standard deviation of the Doppler frequency.

Figure 5 Fig. 5 .
Fig. 5. (a) Mean Doppler frequency as a function of flow velocity angle α.Blue and green curves represent a system with equal Gaussian modes (GG) and a system with Bessel illumination mode and Gaussian detection mode (BG) respectively.Asterisks correspond to values obtained by the numerical simulations shown in panels c and d, and solid lines are calculated using Eq.(27) and Eq.(29).Parameters: NA = 0.3, |v| = 1 mm/s, λ 0 = 780 nm, ∆λ = 120 nm, n = 1.33.(b) Standard deviation of the Doppler frequency.(c) Doppler spectra for different flow velocity angles α, calculated by numerical simulation for a system with equal Gaussian modes.(d) Doppler spectra for a system with a Bessel illumination mode and Gaussian detection mode, NA = 0.3.

Fig. 6 .
Fig. 6.(a) Overview of the algorithm used to obtain functional tomograms containing spatially resolved Doppler frequency spectra and flow velocity components.(b) The interference spectra I(λ , x) are acquired in an over-sampled lateral (x) scan.δ x represents the lateral resolution.(c) A window is applied around each x-position in the lateral scan such that the windowed data represent time samples I w (λ ,t) at that position.(d) Spatially resolved Doppler spectra Ĩ(z, f D ) are obtained after spectral background subtraction, λ to k mapping and 2D FFT.(e) The center of mass µ of the measured Doppler spectrum (blue graph) is calculated and the re-centered data are fitted to Eq. (32) (green graph).The standard deviation σ of the Doppler spectrum is calculated using the fitted curve as weighting function.(f) µ and σ are converted to lateral and axial flow components v t and v z using Eq.(27) and Eq.(29), or their appropriate approximations.

Figure 7 (
Figure 7(a) shows a structural tomogram of the capillary, and Fig. 7(b) and (c) show tomograms of the mean µ[ f D ] and standard deviation σ [ f D ] of the Doppler frequency for a flow rate of F = 0.45 ml/h and angle α = 81 • .As can be seen in 7(a), the capillary is somewhat asymmetric, with an elliptic cross-section.The axial diameter is 200 µm, but the lateral diameter is larger, accounting for a larger cross-section.The intensity gradient over the depth of the flow channel seen in 7(a) is due to three reasons.Firstly, although the extended focus maintains the lateral resolution over 400 µm, it does not have constant intensity over the whole depth, a well-known property of Bessel beams[42].Secondly, the Gaussian detection mode employed by xfOCM

Fig. 7 .
Fig. 7. (a) xfOCM tomogram of the capillary used for flow measurements, dynamic range is 20dB.The channel has an angle α = 81 • and the flow rate is F = 0.45 ml/h.The lateral sampling frequency is 16 times the Nyquist frequency.(b) xfOCM tomogram of the mean Doppler frequencies.(c) xfOCM tomogram of the standard deviation of the Doppler frequencies.10 B-scans were averaged for the images in (b) and (c).

Fig. 8 .Fig. 9 .
Fig. 8. Depth profiles at the center of the capillary (α = 81 • ) measured with xfOCM of (a) the axial flow component and (b) the lateral flow component, for different flow rates.Parabolae are fitted to the measurements, assuming v = 0 at the capillary wall.The maxima of the parabolae are extracted and compared with the expected velocity set by the syringe pump under laminar flow conditions for (c) the axial flow velocity component and (d) the lateral component.Error bars represent ± standard deviation over 10 measurements.