Precision analysis and optimization in phase decorrelation OCT velocimetry.

Quantitative flow velocimetry in Optical Coherence Tomography is used to determine both the axial and lateral flow component at the level of individual voxels. The lateral flow is determined by analyzing the statistical properties of reflected electro-magnetic fields for repeated measurements at (nearly) the same location. The precision or statistical fluctuation of the quantitative velocity estimation depends on the number of repeated measurements and the method to determine quantitative flow velocity. In this paper, both a method to determine quantitative flow velocity and a model for the prediction of the statistical fluctuations of velocity estimations are developed to analyze and optimize the estimation precision for phase-based velocimetry methods. The method and model are validated by phantom measurements in a bulk scattering medium as well as in intralipid solution in a capillary. Based on the model, the number of repeated measurements to achieve a certain velocimetry precision is predicted.


Introduction
Common retinal diseases such as diabetic retinopathy, glaucoma and age-related macular degeneration have shown changes to the vasculature or blood flow during their progression [1][2][3][4][5][6]. Optical coherence tomography (OCT) has become a widely used technique for the studies of disease developments and clinical diagnosis while OCT angiography (OCTA) has been proven to show a particular strong contrast between static tissue and moving particles e.g. blood cells [7]. This makes OCTA already a powerful tool to detect morphological changes in the retinal vasculature and formation of abnormal vessels e.g. neovascularization [8][9][10][11]. Because early symptoms are often connected to reduction of healthy vessels [10] or reduced functionality of existing vessels [5] a quantification of flow is becoming increasingly important. OCT velocimetry is still in its infancy in clinical practice, but it is a prime candidate for further research towards the understanding and early diagnosis of ocular diseases.
The first established technique to determine the flow velocity of particles in OCT was Doppler OCT [12][13][14]. Doppler OCT has the disadvantage that only the velocity component parallel to the propagation direction of the light can be measured. Several attempts have been made in the past years to overcome this limitation, among these: Doppler OCT with multiple directions of illumination/detection [15,16], (complex) speckle decorrelation [17][18][19] and Doppler mean and bandwidth analysis [20,21]. Doppler OCT with multiple directions of illumination requires a significant increase in system complexity [15,16], which can be a limiting factor especially if other functionalities should be incorporated. (Complex) speckle decorrelation and Doppler mean and bandwidth analysis requires in general acquisition of large data ensembles for each voxel and velocity estimation [17][18][19]. Therefore, this approach was demonstrated in situations when samples can be fixated and bulk motion can be Vol. 10 suppressed. The long acquisition times required to acquire the large data sets for these methods limits the applicability for ophthalmology.
An important constraint for the applicability of OCT velocimetry in vivo is the limited acquisition time and thus the number of samples at a tissue voxel that are available to achieve a certain flow estimation precision. An analysis investigating the minimally required number of samples and the optimal method or combination of methods to estimate the flow velocity precision has to our knowledge not been presented. Therefore we set out to provide a theoretical derivation of the probability density function (PDF) of changes in the OCT signal between repeated measurements at (nearly) the same location as a result of motion in individual voxels, and derive a prediction of the precision of flow velocity estimations based on the PDF.
Since scan patterns optimized for flow analysis may require different time delays between repeated measurements, we restricted the analysis to PDFs that only describe the signal change between two individual measurements. The motivation for that is the following: some methods [19][20][21][22] require all samples at (nearly) the same location belonging to one ensemble to be a continuous trace of equidistantly sampled data in time and are optimized for M-scanlike acquisition modes. In many methods for in vivo applications various scan patterns are used in which time intervals between scan repetitions are optimized for sensitivity in a particular velocity range or multiple time intervals are used within the same estimation [8,23]. Our analysis focuses on the latter data acquisition method using only pairs of data for a time interval. As a first demonstration of this analysis we choose the use of Doppler OCT in combination with phase decorrelation as previously introduced [20,21,24].
The manuscript is organized as follows: In section two the theoretical model is described, followed by a numerical simulation based on the theoretical model. An established method by Bouwens et al. [21] is compared to additionally proposed methods. The limit of the Cramér-Rao lower bound is predicted and compared to the achieved precision. In section three the experimental setup and the data processing is described. In section four the experimental measurements and theoretically predicted results are compared. The article will be concluded by a short summary and a discussion about potential future extensions of the theory.

Theory
In order to describe light incident on the sample and the collection of backscattered light in the sample arm of an OCT system, an experimental situation similar to the description by Bouwens et al. [21] is considered. An electrical field propagates from the principal plane of a lens or objective to the scattering volume of a sample and back. The incident field on the volume of scattering particles can be described as superposition of (monochromatic) plane waves with weight ( ) in m in p : the scattering strength of the n-th particle in the first order Born approximation, similar to ref [21]. with n the refractive index. When the light is collected by a lens with focal length f , and the focal point of the lens is located at the origin, the Green's function 2 1 ( ; ) , G out r r p for propagation to the principal plane of the lens can be approximated assuming that the focal length f is much larger than the focal volume from which the backscattered light is collected, i.e.
The detected field has to be weighted by the acceptance mode ( ) out m out p with which light is collected, and 2 r is replaced by the focal length f : In case of a fiber based system the numerical aperture and possible propagation modes determine the accepted modes. For simplicity only point scatterers are considered, Now Eq. (7) 6)) as a function of r , respectively. Next we will determine the change in the phase due to flow. In our analysis we ignore the effect of Brownian motion. Velocities and time differences between repeated scans are considered for which we assume Brownian motion has no significant influence on the signal. Furthermore, we assume bulk motion, a collective motion of all particles in the focal volume. Therefore, flow can be modelled as a translation of all scatterers in the focal volume by a vector δ r . Following the approach of Vakoc et al. [24], the effect of flow is described by a displacement in a second measurement with a bulk transversal and/or axial translation of the detection volume. Two consecutive measurements of the sample result in phasor A Γ and phasor B Γ . The second phasor B Γ can be decomposed into a contribution of the first phasor A Γ and a random phasor C Γ representing the decorrelation. Phasor B Γ can then be expressed as: where A μ and C μ are parameters reflecting the weight of A Γ and C Γ in the composition of B Γ and should adhere to the normalization 2 The weight A μ can be set equal to a parameter α ( A μ α = ), where Vakoc demonstrates this α to be equal to the lateral overlap integral of the detected fields of consecutive measurements. Here we evaluate the volume overlap integral including the axial coherence gate. The field reflected from detection volume A is evaluated as, The field phasor B Γ reflected from detection volume B after motion is defined by the following substitutions: z z z δ → + for axial motion and x x x δ → + for lateral motion. The overlap integral is defined as: Equation (11) cannot be solved analytically for the most general case but we will demonstrate that a solution can be found under the condition that the Rayleigh length R z is much larger than the coherence length c l . In addition to the displacement by flow we also evaluate the effect when the sample is significantly out of focus. For this purpose an additional displacement in z is introduced: the substitution z z z → − Δ is performed in A E . Since the same location but with motion of the scattering particles is described the The numerator of Eq. (11) can then be written as: Now the independence of the overlap integral on the defocus will be demonstrated, i.e., the independence on z Δ . First the effect of z δ has to be evaluated. Due to a rather rapid decorrelation caused by the axial resolution (first line of Eq. (12)) the effect of z δ on ( ) w z and ( ) R z for low numerical apertures will be negligible. The term of the Gouy-phase in the second line and the denominators in the exponents of the third line are then considered constant and independent of z δ . Next the substitution z z z ′ = −Δ can be made in Eq. (12), resulting in: Now we will turn to the integration over z ′ . In most common OCT applications the coherence length c l is much smaller than the Rayleigh-length to be able to resolve many pixels along z with more or less uniform lateral resolution. The exponent determines the range over which the integral along z ′ receives significant contributions. Within this range of order c l , ( ) w z ′ and ( ) R z ′ do not change significantly and are therefore approximated by constants: ( ) w ξ and ( ) R ξ with ξ being the average value of z ′ . This transforms Eq. (13) into: Now only the integrant in the third line of Eq. (14) depends on z ′ . The integrals are evaluated at the center of the coherence gate ( 0) ξ = in analogy to the derivation previously described by Popov et al. [22]. In conclusion, the explicit solution of Eq. (11) is independent of the location of the coherence gate with respect to the beam focus and is given by Equation (15) describes the effect of flow velocity on the overlap integral if we replace the displacement r δ by v τ ⋅ , with v the velocity and τ the time delay between consecutive measurements. Equation (15) is also used in methods for speckle decorrelation for which this factor describes the decay of the autocorrelation function [18,19,22]. It should be emphasized here that Eq. (15) shows that the overlap integral is independent of the distance of the scattering volume to the focal plane, as derived by Popov et al. [22] and previously demonstrated experimentally by Weiss et al. [28]. The remaining calculation of the expected distribution function of phase differences is analogous to Vakoc et al. [24] with one difference: In [24] only lateral motion was considered and investigated. To achieve the PDF of phasor B Γ , the PDF of the sum of the phasors A Γ and C Γ has to be calculated (Eq. (9)), with the prefactor α determined by the overlap integral. The PDFs of the individual phasors are known to follow a Rayleigh and circular Gaussian distribution, respectively [24,29]: The PDF B P for phasor B results from the convolution of A P and C P . It should be pointed out that due to a difference in definition of the fields in the overlap integral (Vakoc et al. defined only incident fields, we define incident fields and subsequent backpropagation) we have The complex part of Eq. (15) can be identified with the mean Doppler shift caused by axial motion. It represents a circular shift of the final PDF of phase differences by a phase 0 2k z δ . Since in this analysis we are interested in the shape of the PDF and not so much its center position, we neglect the phase 0 2k z δ in the further calculation. The convolution of A P and C P leads to [24]: The goal of any estimator is now to estimate the parameter α given a set of measurements from which flow velocities can be calculated using the ratio of displacement and axial/ lateral resolution as given in Eq. (15).

Simulations
In this subsection we simulate multiple repeated phasor measurements for a given overlap integral value α (and therefore a given flow velocity) to create an ensemble of measurements. Randomized ensembles for the same α are created repeatedly to determine the theoretical estimation precision. In all following graphs the uncertainty plotted as standard deviation of the estimated velocities is used as an inverse measure for precision. Thus, high uncertainty indicates poor precision. Equation (19) is used as a starting point for any simulations within this subsection. Further assumptions are: 1. that similar to Srinivasan et al. [18] an isotropic voxel size is used making the broadening of the phase distribution equally sensitive to axial and lateral motion; 2. during prior processing steps the distribution has been shifted circularly to center the mean phase difference at 0 rad. This eliminates the axial Doppler shift 0 2k z δ . Those two assumptions have the same effect as if only lateral motion is simulated.
One method to extract α from phase sensitive measurements was adapted from Bouwens et al. [21]. The distribution of phase differences was fitted with a Gaussian function: where σ is the standard deviation of the Gaussian function, a is a factor to account for nonnormalized data and h is an offset. The standard deviation is then interpreted as ratios  In the following a simulation of measurements based on a finite number of phase difference samples and the effect on the flow precession for different estimation methods is described. For every velocity and therefore for every ratio 0 / r w δ which is subsequently referred to as normalized velocity a number of samples N is drawn from the distribution. All computations are performed in Matlab 2014a (The Mathworks, Inc.). A set A of N randomized uniformly distributed numbers between 0 and 1 are generated (rand() function). The PDF ( , ) P ϕ α is computed and the cumulative probability function ( , ) C ϕ α is derived. The numbers of set A are projected to a set ( ) D ϕ of phase differences according to ( , ) C ϕ α . Three estimation methods are used to estimate the velocity from the simulated phase difference measurements. 1. Gauss. fit: a histogram with 100 bins of phase differences is created from ( ) D ϕ spanning from the minimum to the maximum phase difference, Eq. (20) is fit to it and σ is related to the normalized velocity as discussed above. Of three determined parameters this method uses only one and therefore excludes information which could be used for a reduction of the estimation uncertainty. In order to use the additional information a metric must be found to combine data from the offset together with the data from the Gaussian curve. This requirement and the similarity to phase variance algorithms motivated the next estimation method. 2. second moment: the standard deviation as statistical second moment is computed. As the second moment also takes into account measurements which in the Gauss. fit contribute only to the offset, it is not linearly related to the velocity anymore. This nonlinear relation was accounted for by a look-up-table (LUT) which was used to assign the second moment to a normalized velocity. The LUT was generated in advance by the evaluation of the second moment of ( , ) P ϕ α for a defined range of velocities. Using the second moment is a method to characterize the PDF. To extract all remaining information, all properties of the PDF must be used which motivates the last estimation method. 3. MLE: the normalized velocity was determined by a maximum likelihood estimator (MLE) defined as: where j is the index of the j -th phase difference. For the implementation used here, the likelihood for each phase difference was determined the following way: prior to the estimations, PDFs were created for a range of normalized velocities between 0.01 and 1.8 in 200 steps where each PDF had a resolution with 1024 datapoints in phase space from -π to π. For a phase difference j ϕ the two nearest phase difference datapoints and the respective values of the PDF were found and the likelihood value was determined by linear interpolation between the values of the PDF. This was performed for every j ϕ in an ensemble for all velocities and the one velocity was chosen as estimated velocity which maximized the joint likelihood function. The estimators were determined over many iterations where a new set ( ) D ϕ was created for each iteration. The accuracy and precision (as statistical fluctuation/ uncertainty) of each estimator was quantified by the mean and standard deviation of their results, respectively. Figure 3 shows the results of this simulation for N = 290, 20 different normalized velocities and 5000 iterations.   Fig. 3b due to high overlap with the MLE. The statistical fluctuations are increasing the strongest for the Gaussian fit method. This is shown more quantitatively in Fig. 3c where the standard deviation of the velocity estimation is plotted as a function of the velocities. Here the Gaussian fit method shows a strong nonlinear increase with velocities. Especially for high velocities it is increasing much faster than for the other two estimation methods. To show the differences between the second moment and the MLE, the result of Fig. 3b is plotted on a smaller vertical axis in Fig. 3c. In Fig. 3c the Cramér-Rao lower bound was included and is defined as [30]: .. E the expectation value of any variable inside the brackets. In Eq. (22) the PDF is ( | ) P ϕ α . Equation (22) describes the lowest standard deviation possible due to the sensitivity of ( | ) p x θ to a given parameter (here the normalized velocity). The simulated precision of the flow estimate of the MLE follows exactly the Cramér-Rao lower bound and deviates slightly from the Cramér-Rao lower bound only for the highest velocities plotted. The second moment method has a small offset in comparison to the MLE. The precision of the second moment method is worse than the Gauss. fit method for relative velocities up to 0.25. Hence, it does not perform better than the Gauss. fit method for small velocities (smaller than 0.25) but represents still a significant improvement for the largest part of the velocity range. The MLE clearly outperforms all other methods. Ideally, the statistical error in the determination of flow velocities should increase only linearly. This is the case for normalized velocities between 0 and 1, indicating the velocity range providing the best velocity estimate.

OCT-system
The experimental system shown in Fig. 4 was described in [31] and was modified for this study to be used not only for polarization sensitive data acquisition but also for phasesensitive and therefore Doppler measurements. The light from a 1-μm wavelength swept source (Axsun Technologies, Inc., MA, USA) was split into a lead for the sample arm and reference arm. After being sent through a polarization delay unit (PDU) it was split in the sample arm again by a 20/80 coupler. With respect to the earlier presented system [31] the calibration mirror in one lead of the 20/80 coupler was replaced by a combination of circulator, fiber-bragg-grating (FBG) and photodiode delivering a stable A-line trigger to the acquisition board (Alazar ATS 9360). 10% of the light in the reference arm was split off to a Mach-Zehnder interferometer for k-clocking. After recombination of sample and reference arm the signal was again split into two polarization channels and recorded via balanced detection.

Phantoms
After the signal was guided through the ophthalmic interface it was sent to a model eye shown in Fig. 5. a). The model eye included a plano-convex lens (15mm focal length, LA1540-B, Thorlabs Inc.) and a scattering medium (0.2% Titanium dioxide in cured silicone) with an embedded glass capillary (inner diameter 150 µm, Molex Inc.). The capillary was connected to a syringe pump (NE-4000, New Era Pump Systems Inc.) to create a constant flow rate of a 0.5% aqueous intralipid solution. The fluid was pumped with a flow rate of 0.371 ml/min. Figure 5b shows the structural image of a single B-scan of the phantom. The glass capillary is embedded at a depth of approximately 300 μm.

Data processing
The configuration of the system for polarization sensitivity creates four images in total by depth-multiplexing in the PDU and detection of two orthogonal polarization channels. A more detailed description can be found in [31]. First, a single phase image must be computed out of the individual images. After initial processing steps of fixed pattern removal and dispersion compensation the displacement between the depth-multiplexed images is identified and the images are overlapped. The four fields are averaged for each A-line according to [32]: 12   is used to account for average phase differences between a reference image (here the indexes 1,1 were chosen) and the other three images. Bulk motion compensation has been applied to ( ) out E z similar to previously described in [33].
The calculation of normalized velocities and the standard deviation of their fluctuation has been processed in two ways. Intralipid: adjacent A-lines (0.77 μm apart) were used to calculate phase differences from which velocities are estimated. Each B-scan has been scanned repeatedly to create ensembles for each pixel. The mean Doppler shift was determined by circular averaging and afterwards the phase distribution was shifted circularly to bring the mean phase difference to zero. Under the condition of laminar flow a centrosymmetric flow profile is expected inside the capillary. Therefore pixels belonging to the same velocities are located on circles around the center or the capillary. The inner wall of the capillary was segmented manually and the isotachs are defined by a reduction of the diameter of the wall segmentation centered at the capillary. To avoid effects caused by multiple scattering such as those known as OCT angiography projection artifacts [34] only velocities on the upper half of each circle were used for the analysis of the precision. Static medium: adjacent A-lines (0.77 μm apart) were recorded with a high degree of beam overlap. By comparing A-lines with different displacements, flow can be mimicked because the static medium will look like a homogeneously moving object. For each displacement multiple measurements were collected in ensembles and the velocity was estimated for each ensemble. Afterwards the mean and standard deviation of the velocities were calculated for each displacement. Figure 6 shows the results for measurements of the static medium (a and c) and intralipid (b and d). According to theory the level of the precision of the velocity estimation depends on the number of samples per estimation. Thus, the measurements have been evaluated for a high number (N = 290, a and b) and a low number (N = 33, c and d) of phase differences. The uncertainty as standard deviation was plotted as a function of the experimentally determined normalized velocities. The results from the second moment are not plotted for the measurements of intralipid for better visibility of the other estimators. Overall, a good to excellent agreement was found between theory and experiment. In all cases, the MLE performs significantly better than the Gauss. fit or second moment method, and the second moment method performs better than the Gauss. fit except for small velocities. Figure 6 demonstrates that the MLE provides a velocity precision that even for a number of phase differences as low as 33 is in close agreement with the theoretical prediction. For a larger number of samples and flow in a capillary (Fig. 6b) the experimentally determined flow precision is poorer than predicted by theory, especially at lower velocities. We attribute this deviation to the precision with which the isotachs could be determined in the flow phantom. Close to the capillary wall where the flow velocity is low, the velocity gradient with radius is greatest, introducing the largest variation in velocity with voxel position.

Results
Finally, we present the velocity precision as a function of required samples for both the Gauss. fit and the MLE. This allows selection of the number of samples required to reach a certain velocity precision. Figure 3c, 3d, Fig. 6a-d all show that the velocity estimation uncertainty increases with relative velocity.  Figure 7 shows an additional simulation where we evaluated the velocity estimation uncertainty at a normalized velocity of 0.92 for the Gauss. fit and MLE as a function of the number of samples per estimation in a double logarithmic plot. The uncertainty was simulated from 8 samples to 400 samples per estimation in 40 steps. For uncorrelated data the velocity estimation uncertainty is expected to be inversely proportional to the square root of the number of samples per estimation. Linear functions with slopes of −0.5 were fit to the simulated uncertainties (dashed lines). The velocity estimation uncertainty follows closely the expected behavior, however for a low number of samples a small deviation can be seen. From the shift along the x-axis between the two fits a factor was derived that indicates by how much the number of samples can be reduced by using the MLE instead of the Gauss. fit achieving the same velocity estimation precision. The MLE achieves the same precision as the Gauss. fit with a factor of 8.9 less samples.

Discussion
In the study the estimation precision of predominantly lateral flow velocities has been analyzed for an established method [21] and two proposed estimators. Both estimators have shown lower statistical fluctuations and therefore better precision for a large range of velocities in comparison to the established method. The MLE delivered the best results but already the use of the standard deviation as a statistical second moment represents a significant improvement. The latter is similar to approaches using phase variance for OCTA (e.g [35].) but a LUT as used in this work or calibration curve is necessary for accurate velocity estimations. A more practical demand on measurements, rather than the improvement of the precision, is that a specific precision is desired and ought to be reached with a minimal amount of measurements in order to keep the data acquisition to within clinical acceptable times. For this demand the more interesting aspect is by how much the number of samples can be reduced by choosing the optimal estimation method. This was tested for an exemplary case when the normalized velocity was set to a fixed value and a reduction in samples between the established method (Gauss. fit) and MLE by a factor of 8.9 was found. For in vivo measurements this can be especially helpful where short acquisition times are important to reduce the influence of motion artifacts.
The analysis presented here is based on phase sensitive velocimetry methods. Other methods have been shown, e.g. based on complex speckle decorrelation [19,22]. In the future we plan to extend the theoretical description and simulations to include the amplitude of the OCT signal and use it for a more general analysis that can be applied to methods using the full field of the OCT signal.
Velocities presented in this work are always reported as normalized displacements but the conversion to actual velocity is straight forward if 0 w and the time delay τ over which phase differences are recorded are known. The determined beam overlap 0 / r w δ between neighboring pixels in a static medium was used to calculate the experimental beam radius as 9.59 μm, which was in good agreement to the theoretical diffraction limit of 8.98 μm. determined experimental flow velocity of the intralipid solution in the capillary was 864 mm/s. Assuming a parabolic velocity profile (laminar flow) the maximum velocity given by the set flow rate was 700 mm/s. The feature of the MLE that the estimation uncertainty is approximately linearly proportional to the velocity itself means that the estimation uncertainty is a fixed percentage of the estimated velocity (up to a normalized velocity of 0.9-1.0). The displacement r δ is a function of velocity v and the time interval τ ( ) r v δ τ = which elapses between pairs of measurements from which phase differences are computed. A combination of time delays can be chosen to cover the full range of flow velocities in e.g., the retina, providing an estimation uncertainty which is a fixed percentage of the estimated velocity over the full range.
The optimization of time intervals between repeated scans for velocities expected in the human retina was demonstrated elsewhere [8,23] as well as an additional combination of different time intervals to extend the dynamic range of the velocity estimations. Alternatively, also the beam diameter could be changed.
The presented theory does not include any effects of particles of different sizes or the rotation of non-spherical objects. However, previous studies have shown that methods for velocimetry developed without including those factors can still remain translatable to measurements of blood particles [18,36]. This may still be important in future and require a higher level of numerical simulations but in this work we focused in the theoretical part on an analytical description of the distributions of the phase differences under influence of flow.
In the experiments presented here, the signal-to-noise ratio (SNR) was in the range of 30-35 dB. According to Park et al. [37] a SNR of 30 dB leads to a phase noise standard deviation of 0.022 rad and therefore approximately to a ratio 0 / r w δ of 0.022 which is significantly below all simulated and measured values. Therefore, we treated the effects to the velocity estimations in this work as negligible. However, for in vivo experiments this contribution will become more important as also areas with lower SNR will be present. Effects of low SNR to the determination of the Doppler frequency in phase sensitive OCT and joint spectral and time domain OCT have been investigated previously [38,39]. In such cases, noise from low SNR will change the PDF and therefore change the uncertainty of estimated velocities and the accuracy of the estimation because under such circumstances the measurements of phases contain partially information about decorrelation and partially about the SNR. Similar effects have earlier been treated in the design of an MLE for the Doppler frequency estimation [39].

Conclusion
The PDF was calculated describing the phase fluctuations expected for phasors reflected from a sample with axial and lateral flow. Using this PDF, the statistical limit for the estimation precision of phase based OCT velocimetry methods has been studied and used for the design of an estimator producing results approaching the Cramér-Rao lower bound. It has been compared to established techniques and an MLE based method using the derived PDF provided a significant improvement. The method is designed to be applicable to measurements of pairs of phases which gives it a larger versatility in comparison to methods e.g. restricted to M-scans. In future it is planned to extend this analysis to the full complex field for phasor based techniques.

Funding
The Dutch Technology Foundation STW (grant number 12822), which is part of the Netherlands Organisation for Scientific Research (NWO); the Ministry of Economic Affairs; the Health~Holland framework; the European Union's Horizon 2020 research and innovation program (654148); LaserLaB Europe and Heidelberg Engineering GmbH.