## Abstract

Quantitative velocity estimations in optical coherence tomography requires the estimation of the axial and lateral flow components. Optical coherence tomography measures the depth resolved complex field reflected 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 different 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 field changes between consecutive measurements and use it to predict the uncertainties for amplitude-based, phase-based and complex algorithms. The uncertainty of the flow estimations follows from a statistical analysis and is determined by the number of available measurements and the applied analysis method. The model is verified in phantom measurements and the dynamic range of velocity estimations is investigated. We demonstrate that phase-based and complex (phasor) based lateral flow estimation methods are superior to amplitude-based algorithms.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

In the course of many retinal diseases, such as age-related macular degeneration, diabetic retinopathy and glaucoma [1–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–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–15] and optimize estimation precision and accuracy [16–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–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.

## 2. 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 $\tau$ are defined as ${\Gamma _A}$ and ${\Gamma _B}.$ ${\Gamma _B}$ can be expressed as the sum of a fraction of ${\Gamma _A}$ and a random phasor ${\Gamma _C}$ [14,33]:

where $\alpha$ is the volume overlap integral between measurement A and B. In the case of continuous motion of the volume, $\Gamma $ 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 $\alpha$. 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]:We address now the derivation of a joint PDF for ${\Gamma _A}$ and ${\Gamma _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 $|\alpha |{\Gamma _A}$ and $\sqrt {1 - \alpha {\alpha ^ \ast }} {\Gamma _C}$ are given by,

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 $\varphi = {\varphi _B} - {\varphi _A}$ is introduced and Eq. (9) is integrated over ${\varphi _A}$:

*q*and $p,$ respectively:

We can verify that the here derived PDFs are in agreement with other theories on speckle decorrelation [27–30,38] which have in common that from the autocorrelation of a time trace of the electrical field *E* the overlap integral $\left\langle {{E_{l + 1}}{E_l}^ \ast } \right\rangle = \alpha$ can be calculated, with *l* being the index for the *l*-th measurement. The same expectation value of ${E_{l + 1}}{E_l}^ \ast $ can be calculated using Eq. (12), giving the same result:

*R*of the scattering particles: $|E |= \sqrt R A.$ 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,\varphi ).$ This issue does not exist for $P(q,\varphi )$ because

*q*is insensitive to any scaling factors in the field amplitude. A second reason is that the numerical evaluation of $P(q,\varphi )$ is more stable than $P(p,\varphi )$ 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 $\varphi ,$ respectively:

*β*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)), $\beta$ can be identified with an axial $(\beta = \delta z/\sqrt 2 {l_c})$ or transversal $(\beta = \delta 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–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 $\alpha$ but the PDFs for $P(q,\varphi ),\;P(q)$ and $P(\varphi )$ remain unchanged. We focus in our analysis on the estimation of $\alpha$ and the conversion to $\beta ,$ 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 $\beta .$ As can be seen in the graphs, the PDFs broaden for an increase in $\beta .$ The phase difference PDF in Fig. 2a is a circular symmetrical distribution around average Doppler shift $\theta$ (here the axial motion was zero, so Doppler shift $\theta = 0rad$) over the interval ${\pm} \pi$. 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.$ Fig. 2c, 2d and 2e show the combined phase difference and amplitude ratio distributions for the same values of $\beta$ as in 2a and 2b. All numerical implementations were evaluated in Matlab 2014a (Mathworks Inc.).

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 $\beta ,$ 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 $\beta .$ 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 ${\sigma ^2} = {({NI} )^{ - 1}}$ with *I* the Fisher-Information [40,41].

*P(x)*the PDF of the variable

*x*and the function

*E*[…] in Eq. (18) the expectation value. This leads to the proportionality $\textrm{CRLB} \propto 1/\sqrt 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.

## 3. Experiments

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

#### 3.2 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 $\beta$ they represent and then sorted in 2024/N ensembles with N phasor pairs per ensemble. Thus, the average and standard deviation of $\beta$ 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 $\beta .$ 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 = \sum\nolimits_k {{E_k}{E_{k - 1}}^ \ast } .$ 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 $\beta ,$ three MLEs were implemented according to:

*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 $\theta$ 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 $\beta$ ranging from zero to a maximum ${\beta _{\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. $\beta$ and $\theta$ are estimated in two steps. In the first step, Eq. (19) is evaluated over a coarse grid of 200 samples for $\beta$ (0 to ${\beta _{\max }}$) and $\theta$ (-π to +π), respectively. It was implemented as a search for zero in the first derivative with respect to $\beta$ of the joint likelihood function (the function has only one extremum):

The implementation of the **amplitude-based MLE** is simpler in comparison to the phase-based and phasor-based MLEs because inherent to its definition the amplitude-based MLE is independent of $\varphi$ and $\theta .$ Similar to the phasor-based MLE it is also evaluated in two steps, first with a coarse search for ${\beta _{MLE}}$ and secondly with a finer grid (only two boundaries) for better accuracy.

## 4. Results

Before the MLEs were applied to measurement of intralipid flow in the glass vessel, the linear dependency of the estimation of $\beta$ 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\pi s\kappa /{\rm K}$ was applied with s the shift in pixels, $\kappa$ the wavenumber index and ${\rm 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 $\beta = 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 $\beta = \delta z/\sqrt 2 {l_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*of 13.3 µm. Lateral shifts (6b) are estimated as $\beta = \delta 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.

_{c}The estimation uncertainties of $\beta$ 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 $\beta .$ Thus, a group of ensembles was created from every depth location from every M-scan. The standard deviation of $\beta$ was then plotted as a function of the estimation of $\beta$ 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 $\beta$ and consecutive results were again averaged. This happened by a binning factor of 9 and is shown in Fig. 7(b).

Figure 7(a) and 7(b) show that all three estimation methods achieve the theoretically predicted uncertainty for large ensembles (N = 200). For measurements *in vivo*, an ensemble of 200 measurements requires a significant amount of time. Therefore, we investigated the effect of using a smaller number of measurements in an ensemble. Estimations were also done for N = 50 in 7(c) and 7(d). Figure 7(c) shows the effect of selecting 50 phasor pairs created from 51 consecutive A-lines in an M-scan (Shuffle OFF), or selecting 50 phasor pairs out of consecutive A-lines randomly selected out of the 2025 A-lines in an M-scan (Shuffle ON). The effect of shuffling is shown in 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 $\beta$ 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 $\beta .$ 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 $\beta .$ Those outliers of estimates for $\beta$ increase mean and standard deviation of the whole group of ensembles.

## 5. 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 ($\beta$ 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 ($\beta$ 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 $\beta$ 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, $\beta$ is for transversal motion a function of ${w_0}.$ the flow velocity *v* and the time delay $\tau$ between repeated measurements, $\beta = v \cdot \tau /{w_0}.$ The adjustment of the sensitivity of OCT to different flow velocities has been shown before through the adjustment of $\tau$ [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 $\tau$ 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 $\beta .$ Thus, for a combined range the error would be an approximately fixed fraction of $\beta$ 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.

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

## Funding

Stichting voor de Technische Wetenschappen (12822); Nederlandse Organisatie voor Wetenschappelijk Onderzoek; Ministerie van Economische Zaken; Health~Holland framework; Horizon 2020 Framework Programme (654148); Heidelberg Engineering.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References

**1. **R. Ehrlich, A. Harris, N. S. Kheradiya, D. M. Winston, T. A. Ciulla, and B. Wirostko, “Age-related macular degeneration and the aging eye,” Clin. Interventions Aging **3**, 473–482 (2008). [CrossRef]

**2. **V. Patel, S. Rassam, R. Newsom, J. Wiek, and E. Kohner, “Retinal Blood-Flow in Diabetic-Retinopathy,” Br. Med. J. **305**(6855), 678–683 (1992). [CrossRef]

**3. **M. Emre, S. Orgul, K. Gugleta, and J. Flammer, “Ocular blood flow alteration in glaucoma is related to systemic vascular dysregulation,” Br. J. Ophthalmol. **88**(5), 662–666 (2004). [CrossRef]

**4. **C. L. Chen and R. K. Wang, “Optical coherence tomography based angiography [Invited],” Biomed. Opt. Express **8**(2), 1056–1082 (2017). [CrossRef]

**5. **X. X. Li, W. Wu, H. Zhou, J. J. Deng, M. Y. Zhao, T. W. Qian, C. Yan, X. Xu, and S. Q. Yu, “A quantitative comparison of five optical coherence tomography angiography systems in clinical performance,” Int. J. Ophthalmol. (Engl. Ed.) **11**(11), 1784–1795 (2018). [CrossRef]

**6. **A. Rabiolo, F. Gelormini, R. Sacconi, M. V. Cicinelli, G. Triolo, P. Bettin, K. Nouri-Mahdavi, F. Bandello, and G. Querques, “Comparison of methods to quantify macular and peripapillary vessel density in optical coherence tomography angiography,” PLoS One **13**(10), e0205773 (2018). [CrossRef]

**7. **J. Tokayer, Y. Jia, A. H. Dhalla, and D. Huang, “Blood flow velocity quantification using split-spectrum amplitude-decorrelation angiography with optical coherence tomography,” Biomed. Opt. Express **4**(10), 1909–1924 (2013). [CrossRef]

**8. **S. B. Ploner, E. M. Moult, W. Choi, N. K. Waheed, B. Lee, E. A. Novais, E. D. Cole, B. Potsaid, L. Husvogt, J. Schottenhamml, A. Maier, P. J. Rosenfeld, J. S. Duker, J. Hornegger, and J. G. Fujimoto, “Toward quantitative optical coherence tomography angiography: Visualizing Blood Flow Speeds in Ocular Pathology Using Variable Interscan Time Analysis,” Retina **36**(Suppl 1), S118–S126 (2016). [CrossRef]

**9. **X. J. Wang, T. E. Milner, and J. S. Nelson, “Characterization of Fluid-Flow Velocity by Optical Doppler Tomography,” Opt. Lett. **20**(11), 1337–1339 (1995). [CrossRef]

**10. **Z. P. Chen, T. E. Milner, D. Dave, and J. S. Nelson, “Optical Doppler tomographic imaging of fluid flow velocity in highly scattering media,” Opt. Lett. **22**(1), 64–66 (1997). [CrossRef]

**11. **J. A. Izatt, M. D. Kulkami, S. Yazdanfar, J. K. Barton, and A. J. Welch, “In vivo bidirectional color Doppler flow imaging of picoliter blood volumes using optical coherence tomograghy,” Opt. Lett. **22**(18), 1439–1441 (1997). [CrossRef]

**12. **Y. H. Zhao, Z. P. Chen, C. Saxer, S. H. Xiang, J. F. de Boer, and J. S. Nelson, “Phase-resolved optical coherence tomography and optical Doppler tomography for imaging blood flow in human skin with fast scanning speed and high velocity sensitivity,” Opt. Lett. **25**(2), 114–116 (2000). [CrossRef]

**13. **B. H. Park, M. C. Pierce, B. Cense, S. H. Yun, M. Mujat, G. J. Tearney, B. E. Bouma, and J. F. de Boer, “Real-time fiber-based multi-functional spectral-domain optical coherence tomography at 1.3 mu m,” Opt. Express **13**(11), 3931–3944 (2005). [CrossRef]

**14. **B. J. Vakoc, G. J. Tearney, and B. E. Bouma, “Statistical Properties of Phase-Decorrelation in Phase-Resolved Doppler Optical Coherence Tomography,” IEEE Trans. Med. Imaging **28**(6), 814–821 (2009). [CrossRef]

**15. **J. Walther and E. Koch, “Relation of joint spectral and time domain optical coherence tomography (jSTdOCT) and phase-resolved Doppler OCT,” Opt. Express **22**(19), 23129–23146 (2014). [CrossRef]

**16. **M. Szkulmowski, A. Szkulmowska, T. Bajraszewski, A. Kowalczyk, and M. Wojtkowski, “Flow velocity estimation using joint Spectral and Time domain Optical Coherence Tomography,” Opt. Express **16**(9), 6008–6025 (2008). [CrossRef]

**17. **J. Walther and E. Koch, “Flow Measurement by Lateral Resonant Doppler Optical Coherence Tomography in the Spectral Domain,” Appl. Sci. **7**(4), 382 (2017). [CrossRef]

**18. **A. C. Chan, V. J. Srinivasan, and E. Y. Lam, “Maximum Likelihood Doppler Frequency Estimation Under Decorrelation Noise for Quantifying Flow in Optical Coherence Tomography,” IEEE Trans. Med. Imaging **33**(6), 1313–1323 (2014). [CrossRef]

**19. **Y. K. Tao, A. M. Davis, and J. A. Izatt, “Single-pass volumetric bidirectional blood flow imaging spectral domain optical coherence tomography using a modified Hilbert transform,” Opt. Express **16**(16), 12350–12361 (2008). [CrossRef]

**20. **A. Akif, K. Walek, C. Polucha, and J. Lee, “Doppler OCT clutter rejection using variance minimization and offset extrapolation,” Biomed. Opt. Express **9**(11), 5340–5352 (2018). [CrossRef]

**21. **R. Haindl, W. Trasischker, A. Wartak, B. Baumann, M. Pircher, and C. K. Hitzenberger, “Total retinal blood flow measurement by three beam Doppler optical coherence tomography,” Biomed. Opt. Express **7**(2), 287–301 (2016). [CrossRef]

**22. **A. Wartak, R. Haindl, W. Trasischker, B. Baumann, M. Pircher, and C. K. Hitzenberger, “Active-passive path-length encoded (APPLE) Doppler OCT,” Biomed. Opt. Express **7**(12), 5233–5251 (2016). [CrossRef]

**23. **H. Spahr, C. Pfaffle, P. Koch, H. Sudkamp, G. Huttmann, and D. Hillmann, “Interferometric detection of 3D motion using computational subapertures in optical coherence tomography,” Opt. Express **26**(15), 18803–18816 (2018). [CrossRef]

**24. **L. Ginner, A. Wartak, M. Salas, M. Augustin, M. Niederleithner, L. M. Wurster, and R. A. Leitgeb, “Synthetic subaperture-based angle-independent Doppler flow measurements using single-beam line field optical coherence tomography in vivo,” Opt. Lett. **44**(4), 967–970 (2019). [CrossRef]

**25. **B. Karamata, M. Leutenegger, M. Laubscher, S. Bourquin, T. Lasser, and P. Lambelet, “Multiple scattering in optical coherence tomography. II. Experimental and theoretical investigation of cross talk in wide-field optical coherence tomography,” J. Opt. Soc. Am. A **22**(7), 1380–1388 (2005). [CrossRef]

**26. **N. Uribe-Patarroyo, M. Villiger, and B. E. Bouma, “Quantitative technique for robust and noise-tolerant speed measurements based on speckle decorrelation in optical coherence tomography,” Opt. Express **22**(20), 24411–24429 (2014). [CrossRef]

**27. **N. Uribe-Patarroyo and B. E. Bouma, “Velocity gradients in spatially resolved laser Doppler flowmetry and dynamic light scattering with confocal and coherence gating,” Phys. Rev. E **94**(2), 022604 (2016). [CrossRef]

**28. **V. J. Srinivasan, H. Radhakrishnan, E. H. Lo, E. T. Mandeville, J. Y. Jiang, S. Barry, and A. E. Cable, “OCT methods for capillary velocimetry,” Biomed. Opt. Express **3**(3), 612–629 (2012). [CrossRef]

**29. **I. Popov, A. Weatherbee, and I. A. Vitkin, “Statistical properties of dynamic speckles from flowing Brownian scatterers in the vicinity of the image plane in optical coherence tomography,” Biomed. Opt. Express **8**(4), 2004–2017 (2017). [CrossRef]

**30. **N. Weiss, T. G. van Leeuwen, and J. Kalkman, “Localized measurement of longitudinal and transverse flow velocities in colloidal suspensions using optical coherence tomography,” Phys. Rev. E **88**(4), 042312 (2013). [CrossRef]

**31. **H. W. Ren, K. M. Brecke, Z. H. Ding, Y. H. Zhao, J. S. Nelson, and Z. P. Chen, “Imaging and quantifying transverse flow velocity with the Doppler bandwidth in a phase-resolved functional optical coherence tomography,” Opt. Lett. **27**(6), 409–411 (2002). [CrossRef]

**32. **A. Bouwens, D. Szlag, M. Szkulmowski, T. Bolmont, M. Wojtkowski, and T. Lasser, “Quantitative lateral and axial flow imaging with optical coherence microscopy and tomography,” Opt. Express **21**(15), 17711–17729 (2013). [CrossRef]

**33. **M. G. O. Gräfe, M. Gondre, and J. F. de Boer, “Precision analysis and optimization in phase decorrelation OCT velocimetry,” Biomed. Opt. Express **10**(3), 1297–1314 (2019). [CrossRef]

**34. **B. Braaf, K. A. Vermeer, K. V. Vienola, and J. F. de Boer, “Angiography of the retina and the choroid with phase-resolved OCT using interval-optimized backstitched B-scans,” Opt. Express **20**(18), 20516–20534 (2012). [CrossRef]

**35. **T. Park, S. J. Jang, M. Han, S. Ryu, and W. Y. Oh, “Wide dynamic range high-speed three-dimensional quantitative OCT angiography with a hybrid-beam scan,” Opt. Lett. **43**(10), 2237–2240 (2018). [CrossRef]

**36. **C. Akcay, P. Parrein, and J. P. Rolland, “Estimation of longitudinal resolution in optical coherence imaging,” Appl. Opt. **41**(25), 5256–5262 (2002). [CrossRef]

**37. **J. W. Goodman, * Statistical Optics*, Wiley classics library ed., Wiley classics library (Wiley, 2000), pp. xvii, 550 p.

**38. **J. Lee, W. Wu, J. Y. Jiang, B. Zhu, and D. A. Boas, “Dynamic light scattering optical coherence tomography,” Opt. Express **20**(20), 22262–22277 (2012). [CrossRef]

**39. **J. Tang, S. E. Erdener, B. Li, B. Fu, S. Sakadzic, S. A. Carp, J. Lee, and D. A. Boas, “Shear-induced diffusion of red blood cells measured with dynamic light scattering-optical coherence tomography,” J. Biophotonics **11**(2), e201700070 (2018). [CrossRef]

**40. **H. L. V. Trees, * Detection, Estimation, and Modulation Theory: Radar-Sonar Signal Processing and Gaussian Signals in Noise* (Krieger Publishing Co., Inc., 1992), p. 646.

**41. **R. A. Fisher, “Theory of statistical estimation,” Math. Proc. Cambridge Philos. Soc. **22**(5), 700–725 (1925). [CrossRef]

**42. **O. Nadiarnykh, V. Davidoiu, M. G. O. Gräfe, M. Bosscha, A. C. Moll, and J. F. de Boer, “Phase-based OCT angiography in diagnostic imaging of pediatric retinoblastoma patients: abnormal blood vessels in post-treatment regression patterns,” Biomed. Opt. Express **10**(5), 2213–2226 (2019). [CrossRef]

**43. **B. Braaf, K. A. Vermeer, M. de Groot, K. V. Vienola, and J. F. de Boer, “Fiber-based polarization-sensitive OCT of the human retina with correction of system polarization distortions,” Biomed. Opt. Express **5**(8), 2736–2758 (2014). [CrossRef]

**44. **M. Salas, M. Augustin, F. Felberer, A. Wartak, M. Laslandes, L. Ginner, M. Niederleithner, J. Ensher, M. P. Minneman, R. A. Leitgeb, W. Drexler, X. Levecq, U. Schmidt-Erfurth, and M. Pircher, “Compact akinetic swept source optical coherence tomography angiography at 1060 nm supporting a wide field of view and adaptive optics imaging modes of the posterior eye,” Biomed. Opt. Express **9**(4), 1871–1892 (2018). [CrossRef]

**45. **S. Martinez-Conde, S. L. Macknik, and D. H. Hubel, “The role of fixational eye movements in visual perception,” Nat. Rev. Neurosci. **5**(3), 229–240 (2004). [CrossRef]

**46. **M. Rucci, P. V. McGraw, and R. J. Krauzlis, “Fixational eye movements and perception,” Vision Res. **118**, 1–4 (2016). [CrossRef]

**47. **M. Sugita, R. A. Brown, I. Popov, and A. Vitkin, “K-distribution three-dimensional mapping of biological tissues in optical coherence tomography,” J. Biophotonics **11**(3), e201700055 (2018). [CrossRef]

**48. **P. Cimalla, J. Walther, M. Mittasch, and E. Koch, “Shear flow-induced optical inhomogeneity of blood assessed in vivo and in vitro by spectral domain optical coherence tomography in the 1.3 mu m wavelength range,” J. Biomed. Opt. **16**(11), 116020 (2011). [CrossRef]

**49. **M. T. Bernucci, C. W. Merkle, and V. J. Srinivasan, “Investigation of artifacts in retinal and choroidal OCT angiography with a contrast agent,” Biomed. Opt. Express **9**(3), 1020–1040 (2018). [CrossRef]

**50. **R. F. Spaide, J. G. Fujimoto, and N. K. Waheed, “Image Artifacts in Optical Coherence Tomography Angiography,” Retina **35**(11), 2163–2180 (2015). [CrossRef]

**51. **A. Bouwens, T. Bolmont, D. Szlag, C. Berclaz, and T. Lasser, “Quantitative cerebral blood flow imaging with extended-focus optical coherence microscopy,” Opt. Lett. **39**(1), 37–40 (2014). [CrossRef]