Optical coherence tomography velocimetry based on decorrelation estimation of phasor pair ratios (DEPPAIR)

: Quantitative velocity estimations in optical coherence tomography requires the estimation of the axial and lateral ﬂow components. Optical coherence tomography measures the depth resolved complex ﬁeld reﬂected from a sample. While the axial velocity component can be determined from the Doppler shift or phase shift between a pair of consecutive measurements at the same location, the estimation of the lateral component for in vivo applications is still challenging. One approach to determine lateral velocity is multiple simultaneous measurements at diﬀerent angles. In another approach the lateral component can be retrieved through repeated measurements at (nearly) the same location by an analysis of the decorrelation over time. In this paper we follow the latter approach. We describe a model for the complex ﬁeld changes between consecutive measurements and use it to predict the uncertainties for amplitude-based, phase-based and complex algorithms. The uncertainty of the ﬂow estimations follows from a statistical analysis and is determined by the number of available measurements and the applied analysis method. The model is veriﬁed in phantom measurements and the dynamic range of velocity estimations is investigated. We demonstrate that phase-based and complex (phasor) based lateral ﬂow estimation methods are superior to amplitude-based algorithms.


Introduction
In the course of many retinal diseases, such as age-related macular degeneration, diabetic retinopathy and glaucoma [1][2][3], the blood flow is altered and it is therefore considered to be an important indicator for abnormal developments in the eye. Optical coherence tomography (OCT) is an interferometric imaging technique to measure the coherent backscattering of photons from biological tissue. OCT angiography (OCTA) is a functional extension which has been very successful in recent years for the visualization [4] of blood flow and was even made available as commercial products [5]. However, OCTA only visualizes flow in binary maps. Because early symptoms are often related to small changes in the morphology and functionality of the vessels, quantification becomes increasingly important. Vessel density maps [6] can be used to monitor smaller vessels and capillaries but also the flow velocities and flow rates would be of great interest. Approaches such as SSADA [7] and VISTA [8] have been introduced to extract information about flow rate from modified OCTA measurements and processing.
Numerous approaches have been implemented to retrieve quantitative flow information (primarily velocity) from the OCT signal. The earliest implementation, Doppler OCT [9][10][11], measured the average frequency shift within A-scans. Significant work has been done to improve flow sensitivity by comparing adjacent A-scans [12] describe noise influences [13][14][15] and optimize estimation precision and accuracy [16][17][18][19][20]. However, Doppler-OCT reveals only information about the axial velocity component. In order to convert to absolute flow velocity the flow direction relative to the incident beam must be known, which can be challenging. Work has been reported to resolve the three-dimensional direction by e.g. multi-beam-illumination [21,22] or synthetic apertures in full field [23] and line field OCT [24] but at the cost of a significant increase of the complexity of the device or other limitations such as increased influence of multiply scattered photons [25]. A set of algorithms based on a statistical approach have been developed in which axial and transversal velocity are linked to the correlation between repeated A-line measurements. Examples are SSADA [7], VISTA [8], Intensity speckle decorrelation [26,27], complex decorrelation [28][29][30], Doppler frequency bandwidth analysis [31,32] and phase decorrelation [33]. Those algorithms can be applied using any point scanning OCT device but at the disadvantage that they often need a large number of measurements for each velocity estimation. Depending on the chosen method they can be applied to the amplitude, phase or the complex signal of OCT. The equivalence of the average estimation of some of those algorithms has been addressed in previous reports [15,17] but to our knowledge a quantitative comparison of the retrieved information content has been missing, so far. Additionally, some of those algorithms rely on the use of specific scan patterns, such as M-scans which require all measurements to be part of a trace with equidistant time intervals between individual measurements. This can be challenging when significant bulk motion should be avoided within those traces (rather than between pairs of measurements). It also has been shown that in order to make the dynamic range of the velocity estimations suitable for physiologically relevant velocities, the time interval between measurements and therefore the scan patterns require a certain variability [8,34,35].
In section 2 of this article we provide a derivation of the probability density function (PDF) of the changes of the complex OCT E-field under the influence of motion in the sample, combining the prediction of changes in amplitude, phase and the combined complex field (amplitude and phase) in one theory. This gives us the basis to compare the information for amplitude-based, phase-based and complex signal based algorithms. The Cramér-Rao lower bound (CRLB) predicts the best achievable uncertainty for these analysis bases. This provides us with a tool to evaluate the efficiency of each algorithm to extract information about flow from changes in the E-field. In order to provide best flexibility for the choice of scan patterns we restrict our analysis to the use of an ensemble of pairs of measurements where a pair is acquired sequentially with a time delay. Our analysis requires only a correlation between the measurements in one pair but not between the pairs. In section 3 the experimental system, measurements and data processing is described. Three maximum likelihood estimators (MLEs) are implemented to make the most efficient use of the available data. In section 4 the theory is compared to experimental data and the dynamic range of the estimation is investigated. In sections 5 and 6 the results are discussed and a conclusion is given. The results are an important step towards the understanding and optimization of flow estimations.

Theory
To model the PDF of E-field changes an approach was chosen which was previously introduced by Vakoc et al. [14]. OCT-measurements of the complex E-field are represented as vectors, which we will refer to as here phasors, in a two-dimensional plane where the x-axis represents the real and the y-axis the imaginary part. The phasors of two measurements separated by a time interval τ are defined as Γ A and Γ B . Γ B can be expressed as the sum of a fraction of Γ A and a random phasor Γ C [14,33]: where α is the volume overlap integral between measurement A and B. In the case of continuous motion of the volume, Γ will show a random walk restricted by the amplitude of the reflected field in the experiment and the average step size in this random walk is regulated by α. A movie showing a simulation of the complex E-field backscattered from a collection of moving point scatterers is shown in Fig. 1. For the scope of this paper we chose to use the approximation of bulk flow i.e. all scatterers in the probed volume move with the same velocity. The volume overlap integral is then given by [33]: with δx and δz the lateral and axial displacements of the scatterers, respectively, w 0 the beam radius and l c the coherence length defined as l c = 2 ln 2/π · λ 0 2 /∆λ [36] with λ 0 the central wavelength and ∆λ the full width half maximum bandwidth. In Eq. (2) the phase of α represents the (average) Doppler shift and the absolute value represents the decorrelation by axial and transversal motion. √ 2l c (axially). Right: the E-field which is a coherent summation of the backscattering and collection of all point scatterers. Each step represents a motion of δx/w 0 =0.5 (see Visualization 1) We address now the derivation of a joint PDF for Γ A and Γ C . For fully developed speckle patterns the two phasors in Eq. (1) follow individually a Gaussian distribution in the complex plane [14,37]. We express the real and imaginary part as coordinates in a two-dimensional plane, and the PDFs of |α|Γ A and √ 1 − αα * Γ C are given by, Here, only the decorrelation effect of the absolute value of α is considered. Later the joint PDF for Γ A and Γ C is transformed into the joint PDF for Γ A and Γ B where also the Doppler shift 2k 0 δz of α is incorporated. Because the variables in Eq. (3) are independent of the variables in Eq. (4), the joint PDF can be expressed as the product of the individual PDFs: In order to turn the PDF for fractions of the phasors into a PDF for the phasors with full length, the substitutions x A = x 1 /|α|, y A = y 1 /|α|, x C = ∆x/ √ 1 − αα * and y C = ∆y/ √ 1 − αα * are made in Eq. (5): In practice, Γ C cannot be measured directly but only the pair Γ A and Γ B . Therefore, the joint PDF for the latter phasor pair is derived with the substitutions Equation (7) depends only on the absolute value of the overlap integral α. Next, we incorporate the Doppler shift or phase term of Eq. (2). The effect of a factor e iθ in the complex plane is a pure rotation around the origin by an angle θ. To translate this effect to real valued vectors, which we use to describe the measurements Γ A and Γ B , a rotation in a two-dimensional plane is generated by a rotation matrix: Note that the rotation is only applied to Γ B but not Γ A because Γ A represents the phasor before the phase shift. Replacing x B and y B in Eq. (8) and changing to polar coordinates yields: with A j the radial and ϕ j the angular components of Γ j (j = A, B) representing field amplitude and phase of the measured phasors, respectively. Equation (9) now contains all information about the PDF for a measured pair of phasors. For practical applications not all information is of interest or even meaningful. Therefore, in the following paragraph parametrizations are introduced which help to derive the marginal PDFs that give information about directly accessible values of phasor changes. First, not the actual values of the phase but only the phase change between measurements can be used to determine motion. Thus, the variable for the phase difference ϕ = ϕ B − ϕ A is introduced and Eq. (9) is integrated over ϕ A : To express changes of the amplitude, the coordinate transform p = A A A B , q = A B /A A was chosen. This allows us to express phasor changes as a conjugate product Γ B Γ A * =pe iϕ or a ratio Γ B /Γ A = qe iϕ . Then Eq. (10) becomes: From Eq. (11) the marginal PDFs for P(p, ϕ) and P(q, ϕ) can be derived by integrating over q and p, respectively: with K 0 the modified Bessel function of the second kind. We can verify that the here derived PDFs are in agreement with other theories on speckle decorrelation [27][28][29][30]38] which have in common that from the autocorrelation of a time trace of the electrical field E the overlap integral E l+1 E l * = α can be calculated, with l being the index for the l-th measurement. The same expectation value of E l+1 E l * can be calculated using Eq. (12), giving the same result: For methods to estimate α the PDF P(q, ϕ) is used since it has an easier practical implementation: A A and A B are currently modeled such that they show a Rayleigh distribution with a mean of √ π/2.
In an actual measurement the amplitude of the E-field will scale with the average reflectivity R of the scattering particles: |E| = √ RA. Therefore, the reflectivity will also influence the product p. Either R has to be estimated and the measurements have to be scaled by its inverse, or it must be included in the PDF P(p, ϕ). This issue does not exist for P(q, ϕ) because q is insensitive to any scaling factors in the field amplitude. A second reason is that the numerical evaluation of P(q, ϕ) is more stable than P(p, ϕ) in further calculations due to the Bessel function.
The marginal PDFs for pure phase decorrelation and pure amplitude decorrelation can be derived from Eq. (13) by the integration over q and ϕ, respectively: with φ = ϕ + θ. Equation (15) is identical to the expression derived earlier [14,33]. We can express α in the parameters β, representing lateral and axial motion and θ, representing the Doppler shift: β can be understood as a term summarizing all influences which can lead to a decorrelation of the OCT-measurement. This can for instance be bulk motion, particle flow as depicted in the movie of Fig. 1, or even brownian motion. For bulk motion where all scatterers in the probed volume move with the same velocity along a specific axis (as assumed for Eq. (2)), β can be identified with an axial (β = δz/ √ 2l c ) or transversal (β = δx/w 0 ) velocity normalized on w 0 and l c , respectively. The presented theory can also be used to describe more complex scenarios. It has been shown that diffusion [28][29][30]38], clutter [28,38] (the presence of moving and stationary particles in the same detection volume) and flow gradients [27,39] can be included in the parameter. Also different focusing and light collection geometries [32] can be included, affecting α but the PDFs for P(q, ϕ), P(q) and P(ϕ) remain unchanged. We focus in our analysis on the estimation of α and the conversion to β, which can be identified as relative motion under the above mentioned conditions of bulk flow.
For an illustration of the behavior of phasor ratios, phase differences and amplitude ratios for bulk flow, the PDFs of Eq. (13), (15) and (16) are plotted in Fig. 2 for different values of relative motion β. As can be seen in the graphs, the PDFs broaden for an increase in β. The phase difference PDF in Fig. 2a is a circular symmetrical distribution around average Doppler shift θ (here the axial motion was zero, so Doppler shift θ = 0rad) over the interval ±π. The amplitude ratio PDF in 2b is an asymmetrical distribution between 0 and ∞ with half of the integral below q = 1 and the other half above q = 1.  When the motion becomes stronger, β becomes larger and therefore the phasors decorrelate more. This behavior is expressed as a broadening of the PDFs.
With an exact expression for the PDFs describing the measurements, the calculation of the CRLB [40] is straight forward. The CRLB predicts the smallest achievable uncertainty (which we quantify here as standard deviation) with which a parameter can be estimated from a measurement.
It depends on the sensitivity of the PDF with respect to a change of the parameter, in our case the normalized motion β, but is independent of the estimation method. Therefore we will use the terms phasor-based, phase-based and amplitude-based methods for the CRLBs which belong to the PDFs of Eq. (13), (15) and (16), respectively. The CRLBs for those methods are plotted in Fig. 3a as a function of β. Because the phasor-based methods use the full information available by the OCT measurement, it is expected that the respective CRLB shows the lowest uncertainty. This is also confirmed by the black line in the graph. Very similar behavior is shown by the CRLB of the phase-based methods. In contrast to that, the amplitude-based analysis performs significantly worse and the uncertainty increases nonlinearly, especially for high velocities. The estimation uncertainty decreases with an increase in the number of measurements N (in our case pairs of measurements) which are collected per ensemble to perform each estimation increase. The variance is given as σ 2 = (NI) −1 with I the Fisher-Information [40,41].
with P(x) the PDF of the variable x and the function E [. . . ] in Eq. (18) the expectation value. This leads to the proportionality CRLB ∝ 1/ √ N. This relation can also be reversed in order to calculate a factor describing how many measurements are necessary to reach the same uncertainty with the phase-based and amplitude-based analysis, relative to the phasor based analysis. The phasor-based analysis was chosen as a reference since it resulted in the lowest uncertainty. Fig. 3b shows the square of the ratio of the CRLBs of the phase-based and amplitude-based analysis with respect to the phasor-based analysis, which expresses how many more measurements are needed by these methods as a function of β to reach the same uncertainty as the phasor-based approach. The phase-based analysis performs nearly as good as the phasor-based analysis with ratios between two (small velocities) and one (high velocities). To reach the same uncertainty the amplitude-based analysis requires an order of magnitude more measurements than what is necessary with the phase-or phasor-based analysis. b) The square of the ratio of the CRLB of the phase-based and amplitude-based methods with respect to the phasor-based method, which expresses how many more measurements are needed by these methods as a function of β to reach the same uncertainty as the phasor-based approach. The purely phase-based analysis performs nearly as well as the phasor-based but the purely amplitude-based needs an order of magnitude more measurements especially in the high velocities rang (β ≈ 1) to reach the same uncertainty. The black horizontal line indicates a ratio of one.

Setup
A handheld OCT-system was used as described by Nadiarnykh et al. [42]. Briefly, a swept laser source (Axsun Inc.) with central wavelength at 1060nm and a repetition rate of 100 kHz was used with a fiber based system. The sample arm included a handheld piece which is designed to be held approximately 1cm above the cornea of a human eye with a beam diameter of 1.2mm, while the person is facing upwards. Instead of an actual eye, a model eye was used consisting of a lens (AC-127-019-C, Thorlabs Inc.) and a retina phantom as it was previously presented in [33]. A schematic drawing of the measurement configuration and a structural B-scan of the phantom is shown in Fig. 4. The phantom consisted of a bulk scattering material (TiO 2 particles in cured silicone) with an embedded glass vessel (inner diameter 150µm). A 0.5% aqueous intralipid solution was pumped by a syringe pump (NE-4000, New Era Pump Systems Inc.) through the vessel to create flow of a scattering fluid. Flow was measured in two different configurations. In the first configuration, the beam was translated over the static scattering medium, mimicking bulk flow. In the second configuration, flow was measured in the vessel. Data acquisition was performed in M-scan mode, for each lateral location 2025 A-lines were recorded, before laterally displacing the beam by 1.1 µm for the next M-scan. 2000 M-scans were acquired to create a B-scan. This generates a field of view of 2.2 mm and a acquisition time of 40.5 s per B-scan. After a phantom measurement, a noise measurement was performed for fixed pattern background removal by blocking the sample arm. The signal-to-noise ratio (SNR) in the phantom was between 20 and 30dB. Fig. 4. a) schematic drawing of the setup. The scan module of an OCT system emitting a collimated beam was mounted above a lens and retina phantom. The phantom consist of a glass vessel embedded in a scattering medium. The capillary was connected to a syringe pump which pumped an aqueous intralipid with a constant flow rate. b) Structural B-scan of the phantom with indicated location of capillary and area of static scattering medium.

Processing
The first A-line processing steps were background subtraction for fixed pattern removal, spectral shaping with a tapered cosine window and Fourier transform. Secondly, individual measurements in the M-scans were combined into phasor pairs. For the application of the MLEs to experimental data, the phasor pairs were first grouped by the β they represent and then sorted in 2024/N ensembles with N phasor pairs per ensemble. Thus, the average and standard deviation of β could be calculated for each group of ensembles. For the static medium, the phasor pairs were created from two A-lines at the same depth location separated by the same step size (a multiple of 1.1 µm) to represent a group of ensembles with the same β. For measurements in intralipid each M-scan is the measurement of 2025 A-lines at the same location resulting in 2024 phasor pairs for each depth location which represents one group of ensembles.
Each phasor pair was corrected for bulk motion in the following way: the conjugate product was summed along each A-line in the bulk phantom: E = k E k E k−1 * . The angle of this sum was extracted and used as bulk phase shift for compensation. The glass vessel was embedded below a layer of bulk scattering medium. The bulk phase shift which was extracted in an A-line from the medium above the intralipid solution, was also applied to measurements of the intralipid. For the estimation of β, three MLEs were implemented according to: where j={1,2,. . . , N} is the index of the j-th measurement. Equation (19) is the definition for the phasor-based MLE using the PDF of Eq. (13). For the phase-based and amplitude based MLEs it must be replaced by the PDFs in Eq. (15) and (16), respectively. In the following paragraph we describe the implementation of the individual MLEs. The phase-based MLE was described in [33]. In brief: Before applying the MLE, the average axial Doppler shift θ in the phasor pair distribution was removed by circular averaging of the phase differences and afterwards the phase distribution was shifted circularly to bring the mean phase difference to zero. PDFs were created for β ranging from zero to a maximum β max of 2.5 in 200 steps. These PDFs were used to compute the joint likelihood function. In this implementation we added a spline interpolation of the joint likelihood function for a more accurate determination of the maximum likelihood. For the phasor-based MLE the estimation of the average axial Doppler shift is included in the MLE. β and θ are estimated in two steps. In the first step, Eq. (19) is evaluated over a coarse grid of 200 samples for β (0 to β max ) and θ (-π to +π), respectively. It was implemented as a search for zero in the first derivative with respect to β of the joint likelihood function (the function has only one extremum): The numerical evaluation of Eq. (20) on a grid provides pairs of θ and β for which Eq. (20) is close to zeros. These pairs trace out a path in a two dimensional (θ, β) space (see Fig. 5).
Only for the smallest β on this path a unique value for θ can be found. Finding this unique combination is equivalent to finding parameters which maximize the likelihood function and therefore yield β MLE and θ MLE . In a second step the accuracy of the estimation is improved: the four combinations of θ and β around the zero-crossing for β MLE and θ MLE are used as boundaries to compute a second grid with 200 steps for each parameter. The implementation of the amplitude-based MLE is simpler in comparison to the phasebased and phasor-based MLEs because inherent to its definition the amplitude-based MLE is independent of ϕ and θ. Similar to the phasor-based MLE it is also evaluated in two steps, first with a coarse search for β MLE and secondly with a finer grid (only two boundaries) for better accuracy.

Results
Before the MLEs were applied to measurement of intralipid flow in the glass vessel, the linear dependency of the estimation of β on axial and lateral motion was verified and the dynamic range was investigated. This was done by controlled axial and lateral bulk motion of the static scattering medium. To mimic axial motion, a B-scan was created by selecting one A-scan out of every M-scan. Axial motion was generated digitally -a copy of the complex valued B-scan was Fourier transformed to spectral domain, a phase ramp 2πsκ/K was applied with s the shift in pixels, κ the wavenumber index and K the number of samples per A-line. The shifted image was back-transformed to depth domain and phasor pairs were assembled from corresponding pixels in the original and shifted copy of the B-scan. For lateral motion phasor pairs were generated from adjacent A-lines, mimicking a lateral translation of the sample. Due to a high degree of beam overlap, motion in fractions of the lateral resolution could be mimicked. The results for axial and lateral motion estimation are shown in Fig. 6 for the phasor-based MLE for different numbers (N) of phasor pairs per ensemble. The fit of a line shows the linear regime in both graphs (a and b). The upper limit of the linear regime, which defines the dynamic range depends on N. In Fig. 6(a) where the axial motion was mimicked, the linear regime extends on the lower end to zero. In Fig. 6(b) where the lateral motion was mimicked using adjacent A-lines, an offset can be noticed. Here, a line cannot be fitted which crosses the origin. Also close to the origin the measurements deviate from a linear slope. We attribute this deviation to the influence of shot noise and beam position jitter originating from the galvanometer scanners. The lowest measurable β = 0.098 represents the lower bound of the dynamic range and was obtained with a stationary beam (except for galvanometer scan jitter). The slopes in Fig. 6 can be used to determine the coherence length and beam radius. According to Eq. (2) and (17) axial shifts (6a) are estimated as β = δz/ √ 2l c . Assuming a Gaussian spectrum of the source the coherence length l c was 14.6 µm. A fit to the slope in Fig. 6(a) resulted in a coherence length l c of 13.3 µm. Lateral shifts (6b) are estimated as β = δx/w 0 . The theoretically calculated diffraction limited beam radius in the sample is 11 µm, given by the beam radius at the last lens before the phantom and the focal length of the lens. The actual beam radius derived from the slope in 6b is 15 µm.
The estimation uncertainties of β and the agreement with the CRLBs was verified by a measurement with intralipid in the glass vessel. The solution was pumped with a flow rate of 435µl/min. In every M-scan only measurements inside the vessel were used and every location along the depth inside the vessel represents a different β. Thus, a group of ensembles was created from every depth location from every M-scan. The standard deviation of β was then plotted as a function of the estimation of β averaged over the ensembles. To reduce effects of multiple scattering only the upper half of the vessel was analyzed. The results are shown in Fig. 7(a) for an ensemble size of N = 200. The estimations of the measurements are displayed as dots while the CRLBs for the corresponding methods (phasor-based, phase-based, amplitude-based) are displayed as solid lines in the same color as the corresponding MLEs. Figure 7(a) shows that the measurements follow approximately the predicted CRLBs. To increase the number of ensembles from 10 to 90, binning was applied. For the binning the results were sorted according to the average estimations of β and consecutive results were again averaged. This happened by a binning factor of 9 and is shown in Fig. 7(b).   Fig. 7(c) as an example for the phasor-based estimations. Red dots represent estimations from phasor pairs as consecutive neighbors (no shuffling), black dots represent estimations from shuffled phasor pairs. Figure 7(c) demonstrates that the estimations from shuffled phasor pairs achieve lower uncertainty. This example demonstrates that phasor pairs constructed from 51 consecutive M-scans are not completely statistically independent. For all other sub-figures in Fig. 7 we used Shuffle ON.
While the phasor-and phase-based MLE achieve the uncertainties predicted by the CRLBs, large artifacts were found for the amplitude-based MLE in 7d (delineated in black). It was found that those artifacts have a higher chance of appearance for large β and become reduced if N is increased. They exhibit a patch-like organization, indicated by the black ellipse. Each plotted point in the graph represents one group of estimations. In the first patch (bottom left in the ellipse) every group of ensembles contained one ensemble which had a constellation of measurements that resulted in a joint likelihood function for an unexpectedly high β. In the second patch two outliers can be found in the group of ensembles, in the third patch three, etc..
If N is increased, more measurements are included in the ensemble, the expected PDF seems to be better resembled by the measurements and the joint likelihood function peaks again at a more reasonable β. Those outliers of estimates for β increase mean and standard deviation of the whole group of ensembles.

Discussion
In this study a new theory was developed to compare quantitatively the precision of flow estimation by analyzing the information content of the estimation of decorrelation and related motion by three categories of methods: (complex) phasor-based, phase-based and amplitude-based analysis. The PDF's describing the decorrelation between two measurements were derived for these three methods. Based on the PDF's, the CRLB was calculated, and experimental data was compared to the CRLB prediction. The most striking result is that the information on decorrelation is dominated by information from the phase over the amplitude of the OCT-signal. This is especially the case for the upper range of velocities / motion (β close to 1). Because complex methods combine both information from phase and amplitude they have the best theoretically achievable precision. A maximum likelihood estimator was implemented for each of the above mentioned categories and was applied to measurements of a flowing intralipid solution. An excellent agreement of the experimental measurements and theory was found for the complex phasor-based and phase-based methods. For the amplitude-based MLE we found artifacts for the combination of strong decorrelation (β close to 1) and small numbers of samples per ensemble. The dynamic range of the flow estimation was investigated and it was shown that at the upper end it depended on the number of samples per ensemble. For the numbers used here (N = 25 to N = 200) an upper range of β of approximately 0.95 to 1.15 was found. The lower bound was independent of the number of samples and only determined by noise sources and was in the measurements here at approximately 0.1.
The dynamic range of the velocity calculation is an important aspect for applications to measure e.g. blood flow. As presented here, β is for transversal motion a function of w 0 . the flow velocity v and the time delay τ between repeated measurements, β = v · τ/w 0 . The adjustment of the sensitivity of OCT to different flow velocities has been shown before through the adjustment of τ [8,35,43] This can be realized by the adjustment of scan patterns. A typical value for w 0 is 10µm. Under this assumption Table 1 shows three examples how τ can be designed to cover different flow regimes. This can be realized by the adjustment of scanning patterns. Alternatively, also scanning with multiple beam diameters might be possible as it has been reported before [44].  Figure 8 illustrates the combination of dynamic ranges as mentioned in Table 1 with the relative error (standard deviation divided by expectation value) for a number of N = 25 phasor pairs. One interesting aspect for the combination of multiple time intervals to cover a larger range of velocity estimations is that for the phase-based and phasor-based methods the statistical error is nearly linearly dependent on the expectation value of β. Thus, for a combined range the error would be an approximately fixed fraction of β over three decades of flow velocities. The amplitude-based analysis on the other hand results in strong discontinuities in the regions where two ranges are merged. The limiting noise sources for the lower end of the dynamic range are shot noise and motion artifacts. Shot noise can be related to the SNR of a measurement and has been shown to broaden the distributions of Doppler phase shifts [13] and Doppler frequencies [15] and the implementation of both shot noise and decorrelation noise for optimal average Doppler frequency estimation has been demonstrated [18]. In the future, an incorporation of the influence of shot noise into the model described here and estimation of decorrelation would be beneficial because it is expected that it can bias the estimation and influence the uncertainty. A modeling of bulk motion artifacts is more complicated because it can come from different sources e.g. positioning error (jitter in scanner position) and sample motion and will highly depend on the technical implementation and conditions of the measurement. OCT is often applied in ophthalmology where it is used for minimally invasive monitoring and diagnosis of diseases and where sample fixation is very limited. Eye motion as micro-saccades, drift and tremor [45,46] can significantly influence the repeatability of the measurement of the same location. Limiting the data analysis to ensembles of measurements where only pairs of measurements are correlated (as implemented here) rather than the requirement of traces of measurement might help to overcome this limitation.
A fundamental assumption for the theory presented here is that fully developed speckle are measured. A consequence is that measurements follow a random, two-dimensional Gaussian distribution in the complex plane and therefore the amplitude follows a Rayleigh distribution. It has been shown in the past [47] that this is not generally fulfilled in biological tissue. Then the amplitude distribution can be approximated by a K-distribution and different cell and tissue types can be identified by the fitting parameters of this distribution. This changes the speckle pattern and it can be expected that the distributions of phases and phase changes may be altered by this, as well. Since for this study only small scattering particles were used (TiO 2 , intralipid) this effect did not play a role but in future an investigation of this effect to blood flow measurements might be important. The backscattering of blood cells is different from TiO 2 and intralipid and depends on the orientation of the cells [48,49]. Strong forward scattering can also influence the speckle pattern of an A-line profile. Light which is scattered in tissue below a blood vessel can carry amplitude and phase changes due to forward scattering in the blood vessel above.. In this theory the effects of multiple scattering was not included. Those can be important in the explanation of shadowing artifacts known in OCT angiography [50] and the consideration might be important for accurate velocity and flow rate quantification when multiple vessels are within one A-line profile. However, previous reports showed the applicability of similar conditions to blood flow measurements [7,28,39,51] which might indicate that estimations are still valid or only small changes are necessary.

Conclusion
A PDF was derived to describe the changes of the complex OCT E-field between individual measurements under a decorrelation influence such as axial and lateral flow. From this also the PDFs for pure phase and pure amplitude changes were derived. Using those PDFs the CRLBs for the individual categories of analysis methods were derived, predicting the uncertainty of the estimations of decorrelation under the condition of a given number of measurements. MLEs have been implemented for each methods category and were used to verify the theoretical predictions of the CRLBs. The most interesting result is that the estimation uncertainty is almost completely determined by the phase of the measurement while the respective MLE significantly outperforms the amplitude-based MLE. The method requires correlation between measurements in ensembles of pairs but does not require correlation between the individual pairs which provides high flexibility in the choice of scan patterns.