Blood flow rate estimation in optic disc capillaries and vessels using Doppler optical coherence tomography with 3D fast phase unwrapping

The retinal volumetric flow rate contains useful information not only for ophthalmology but also for the diagnosis of common civilization diseases such as diabetes, Alzheimer’s disease, or cerebrovascular diseases. Non-invasive optical methods for quantitative flow assessment, such as Doppler optical coherence tomography (OCT), have certain limitations. One is the phase wrapping that makes simultaneous calculations of the flow in all human retinal vessels impossible due to a very large span of flow velocities. We demonstrate that three-dimensional Doppler OCT combined with three-dimensional four Fourier transform fast phase unwrapping (3D 4FT FPU) allows for the calculation of the volumetric blood flow rate in real-time by the implementation of the algorithms in a graphics processing unit (GPU). The additive character of the flow at the furcations is proven using a microfluidic device with controlled flow rates as well as in the retinal veins bifurcations imaged in the optic disc area of five healthy volunteers. We show values of blood flow rates calculated for retinal capillaries and vessels with diameters in the range of 12–150 μm. The potential of quantitative measurement of retinal blood flow volume includes noninvasive detection of carotid artery stenosis or occlusion, measuring vascular reactivity and evaluation of vessel wall stiffness. © 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
Volumetric blood flow (blood flow rate) in the human eye is one of the parameters that characterize blood circulation. It has a wide range of potential applications in the diagnosis of systemic diseases, e.g., diabetes, cerebrovascular disease, hypertensive heart disease, and Alzheimer's disease, with possible links to ocular conditions such as diabetic retinopathy, vessel occlusion, and hypertensive retinopathy [1][2][3][4][5][6]. Ocular blood flow can also be used to analyze brain perfusion, since the flow to the eye provides direct insight into the internal carotid artery flow dynamics, which provides the majority of the blood flow through the brain. The dynamics of the flow in ocular vessels can be translated into non-invasive flow-based measurement of blood pressure in the cerebrovascular system using a modern computational fluid dynamics approach to signal modeling [7]. Quantitative measurement of blood flow in the brain is currently either expensive, e.g., in the case of magnetic resonance imaging, or invasive, e.g., in the case of computed tomography perfusion. Consequently, this parameter of brain function is rarely considered in routine evaluations. Therefore, fast and noninvasive optical methods for measuring the ocular blood flow rate combined with creative approaches to flow signal processing are required.
Optical coherence tomography (OCT) is a commonly used optical method for 3D structural imaging of the human retina. It has the advantage of providing access to the phase of scattered light, allowing for the detection of small displacements or flows. Therefore, it can be used to provide information about retinal circulatory systems. For example, OCT angiography (OCTA) methods generate images of the blood flow indicating the locations of vessels within the tissues [8][9][10][11][12] and Doppler OCT methods and other OCT velocimetry techniques provide information about blood flow velocity or rate [13,14].
Blood flow rate is defined as the volume of blood flowing through the vessel lumen in a given period of time. There are two approaches to calculating blood flow rates from Doppler OCT data. One requires the knowledge of the mean flow velocity along the vessel and the area of the lumen perpendicular to the blood flow. Multiplication of the two gives the mean flow rate. To determine the volumectric flow the velocity along the vessel has to be measured for each data point across the vessel lumen to account for the change in the velocity from the center to the vessel wall. The flow rate is the integral of the products of velocities and the pixel areas taken over the vessel lumen. The difficulty in this approach lies in the assessment of the vessel orientation in the three-dimensional data set. Several methods have been developed to measure this angle, including segmentation of the vessels in 3D data sets [15], use of circular scan-patterns to determine the local orientation of the vessels [16,17], and two-or three-beam scanning for improved angle measurement accuracy and simultaneous calculation of multiple components of velocity vector [18][19][20]. The second approach is based on the observation that the flow rate can be calculated from a velocity vector projection in any direction as long as the area of the vessel cross section perpendicular to that direction is also known [21]. In Doppler OCT, the axial projection (Z-component) of the blood flow velocity vector is measured. Multiplication of the Z-component of the velocity and the XY-area of the pixel at which the velocity is measured gives the local flow rate. If the vessel lumen is fully visible in the XY-cross section (C-scan), then the sum of the products of the Z-flow velocity and the XY-pixel area calculated at all the XY-positions in the vessel cross section gives the blood flow rate. In this approach, knowledge of the vessel orientations in the 3D data set is not required. However, segmentation of the vessels in the C-scans is necessary if the flow rate needs to be calculated for each individual vessel. Alternatively, a net flow rate through an entire C-scan of interest can be calculated by summing the products of all the local Z-velocities and XY-pixel areas without the need for vessel segmentation. This approach could be of interest in the analysis of blood flow in the optic nerve head.
The blood flow rate can be calculated from Doppler OCT data provided that the axial flow velocity component can be measured. This may be problematic if the blood flow is too fast compared to the OCT imaging rate. Whenever the axial flow velocity component exceeds the maximum velocity measurable by a given OCT system, the phase change of the complex signal representing the flow exceeds the interval (-π, π) and phase-wrapping occurs. The phase wrapping renders characteristic images of flow velocity in the form of concentric rings of opposite flow directions within vessel lumens [22].
Phase unwrapping methods aim to detect the phase wrapping and correct the data such that phase discontinuities (i.e., phase wraps) are removed without the introduction of phase errors. Typically, the correction requires addition or subtraction of multiples of 2π when the phase value in a group of neighboring voxels exceeds the -π or +π limit. Many phase unwrapping methods have been introduced to Doppler OCT [23][24][25][26][27][28]. Recently, we developed a robust and computationally efficient phase unwrapping method, namely four Fourier transform fast phase-unwrapping (4FT FPU) [29], and showed that it is specifically suited for multidimensional Doppler OCT data. Similar method was also presented for 2-dimensional data acquired using stroboscopic interferometry [30].
The focus of our previous article was on the development of algorithms suitable for Fourierdomain OCT imaging [29]. Capability to unwrap the phase was tested in numerical simulation of flow and in a flow phantom for 2-dimensional and 3-dimensional data sets was presented. We have introduced a metric to quantify the phase unwrapping errors and tested the performance of phase unwrapping methods in data with increasing noise level. This analysis enabled us to predict at which phase noise levels the algorithms can be used without introducing errors resulting in incorrect axial flow velocity values. Finally, we have used the phase unwrapping algorithms in Doppler OCT data acquired in the human retina to demonstrate their applicability in calculation of axial blood flow velocity. We have identified the three-dimensional version of the fast phase unwrapping algorithm as better suited for biomedical applications due to its higher immunity to noise than the two-dimensional FPU method.
Although the previous work demonstrated the applicability of the 3D 4FT FPU method to calculate the axial flow velocity in the human retinal vessels, the question still remains whether calculations of volumetric blood flow give reliable results. That is, if they provide correct values of the blood flow on which medical doctors could rely while diagnosing diseases and to what extent they can rely on Doppler OCT measurements. Although it may seem that there is a little difference between calculating the axial flow velocity and volumetric flow rate, there is a considerable difference. Axial component of the flow velocity vector is highly susceptible to variations caused by changing angular orientations of vessels. For example, they may be a result of the anatomy of the vascular system, or alignment of the eye in the OCT device. Because the axial flow velocity does not necessarily reflect only changes in the blood circulation, it is not the best metric for the use in ophthalmic diagnostics. On the other hand, volumetric flow rate, when possible to calculate, is not dependent on the angular orientations of vessels. Although it is not free from possible errors, it is a better choice of metric for ophthalmic applications.
In this article we focus on two questions. First, does the Doppler OCT with the 3D 4FT phase unwrapping method give reasonable values of volumetric blood flow in retinal vessels? Second, can the algorithms be implemented to provide results in real-time or at least within a few minutes after acquiring the data?
Answering the first question isn't an easy task since there is no gold standard to which we could relate our results. That is, there is no scientific consensus as to what are the reference values of blood flow in retinal vessels. Moreover, in the diagnostic practice, medical doctors rarely look for the "ground truths". They relate the results from a specific diagnostic technique or even a specific type of a device to the established norms and pathology classifications. However, only a few papers exist reporting volumetric blood flow values measured with OCT techniques in a limited number of cases [15] [31][32][33][34][35]. Therefore, to answer the first question we have performed two tests, which do not require comparisons with established norms, yet show the reliability of the calculated volumetric flow rate values. In the first test, we check if the additive property of the volumetric flow rate at furcations is fulfilled in the phase-unwrapped Doppler OCT data obtained in a flow phantom (a microfluidic device) and in human retinal vessels (five volunteers). The experiment in the flow phantom also enables us to check if the measured flow values agree with the flow set in the microfluidic pump, thus indicating if phase-unwrapped Doppler OCT generates any offsets in the measured values. In the second test, we measure the blood flow in the human retinal vessels (three volunteers) in a large number of vessels (600) with lumen diameters ranging from 12 to 150 µm. We expect that there should be a well-defined dependence between the vascular lumens diameters and the volumetric blood flow rates [36]. Lack of such dependence would indicate erroneous results of blood flow calculation in the phase-unwrapped Doppler OCT data.
The answer to the second question is important from the perspective of practical clinical applications of the phase-unwrapped Doppler OCT imaging. Efficient and wide-spread use of these methods depends not only on providing reliable results but also on the availability of the Doppler OCT images and blood flow rate measurement results immediately after the data is acquired in the eyes of visiting patients. To demonstrate that real-time phase-unwrapped Doppler OCT imaging is possible, we have implemented the entire data processing pipeline from basic data processing, through implementation of Doppler OCT methods, to phase unwrapping using graphic processing units (GPU) technology.

Materials, subjects, experimental setup, and Doppler OCT imaging methods
Doppler OCT imaging was performed in a flow phantom and in the retinas of five volunteers. The flow phantom was a microfluidic device. A 300 × 100 µm rectangular microchannel was trifurcated into three 100 × 100 µm rectangular channels in a polydimethylsiloxane (PDMS) substrate mixed with titanium dioxide (TiO 2 ) powder, which acted as a light scattering medium. A protein solution consisting of 0.1% cow milk was pumped through the microfluidic device with a syringe pump (neMESYS 290N, CETONI GmbH).
Imaging of the flow phantom was performed with a research-grade spectral-domain optical coherence microscopy (OCM) setup. Light with a center wavelength of 800 nm was emitted by a supercontinuum source (SuperK EXTREME/FIANIUM, NKT). The axial imaging resolution was ∼3 µm. The lateral imaging resolution was 8 µm. The A-scan rate was set to 70 kHz by adjusting the repetition time of the camera used in the OCT spectrometer (Basler Sprint spL2048-140km, Basler AG). The imaging protocol (a) is given in Table 1. Imaging of the eyes of healthy volunteers was performed in compliance with the Declaration of Helsinki [37] and regulations of the Bioethics Committee of the Nicolaus Copernicus University, Torun, Poland. The fundus of five subjects (age range, 26-40 years) was imaged in the optic nerve head area.
The imaging of volunteers was performed with a commercial spectral-domain OCT device (REVO NX, Optopol Technology, Poland). The light source was a superluminescent diode (center wavelength, 830 nm). The axial imaging resolution in tissue was ∼5 µm and the lateral imaging resolution was 12 µm. The A-scan rate was 110 kHz. The imaging protocols (b-c) are given in Table 1.
Doppler OCT imaging in the flow phantom was performed with phase-resolved Doppler OCT method [13] and imaging in the human eyes performed using the joint spectral-and time-domain OCT (STdOCT) method [13]. In Doppler OCT imaging, a series of A-scans is collected from the same spatial location in the object and, in practice, with the change in the beam position between two consecutive OCT spectra acquisitions no larger than half of the beam spot diameter [38,39]. To calculate the blood flow in the human retina, we used groups of four A-scans. The imaging protocols (b-c) are given in Table 1.

Volumetric flow rate calculations
Light back-scattered at moving scaterrers, e.g., flowing blood cells, is received by the detector with a frequency shifted by ω D owing to the Doppler effect. For each spatial location in the sample r, this shift is proportional to the velocity component parallel to the light beam, ω D (r) = 2knv z (r), where k is the central wavenumber, n is the index of refraction of the moving medium, and v z is the axial component of the velocity vector v. In phase-resolved Doppler OCT (PR-Doppler OCT) imaging, the approximation of the first derivative of the Doppler frequency, ϕ D (r, ∆t) = ∆ω D (r)/ ∆t, is obtained using the measured phase differences between pairs of A-scans, A(r, t) = |A(r, t)| exp(jknv z (r)t), delayed by ∆t [40]: In the STdOCT method, the Doppler frequency ω D is obtained directly by Fourier-transforming groups of A-scans along the time coordinate, which allows for direct use of multiple A-scans to increase the signal-to-noise ratio (SNR) and improve the image quality ( Fig. 1) [13]. Knowledge of the Z-component of the flow velocity (v z ) and the vessel cross section S perpendicular to the light beam allows for the calculation of the volumetric flow rate Q according to [21] where ψ(r) is the unwrapped phase difference ϕ D (r, ∆t), s is the pixel area in the cross section, and ROI is the region of interest for which the flow rate is calculated. In practice, ROI contains the cross section of the vessel lumen perpendicular to the light beam as shown in Fig. 3(a). To obtain ψ(r), a phase unwrapping algorithm has to be implemented.

Phase unwrapping
The phase ϕ D (r, ∆t) measured using the PR-Doppler OCT or STdOCT method is always in the range of −π to π. Phase wrapping occurs when the flow velocity exceeds this limit [22]: Measurement of velocities exceeding the v z max value requires implementation of numerical phase unwrapping methods. We used our recently developed [29] three-dimensional (3D) fast phase unwrapping (FPU) calculated using four complex Fourier transformations (4FT) (i.e. 3D 4FT FPU). The FPU method calculates the phase correction ψ est , in the phase differences distribution ϕ, which is subsequently used to calculate the unwrapped phase distribution ψ(r): The Laplace operators ∇ −2 and ∇ 2 are implemented using Fourier transformations [29,30,41].

Implementations of algorithms
We have implemented a processing pipeline consisting of the transformation of OCT spectra into complex A-scans [42,43] followed by one of the two Doppler OCT methods (PR-Doppler OCT [22] or STdOCT [13]), with one of the two versions of the 4FT FPU algorithms [29]. We used the STdOCT method with a window of four A-scans and 128-point Fourier transforms. The 2D and 3D 4FT FPU algorithms were implemented separately and called depending on the dimensionality of the processed data (2D/3D).
The software was developed under the Microsoft Windows 10 operating system. The graphical user interface (GUI) was prepared in LabVIEW 2018. All the processing cores (STdOCT, PR-Doppler OCT, 4FT FPU 2D, and 4FT FPU 3D) were written in C using Microsoft Visual Studio 2017 and compiled into the DLL library. We used a graphics processing unit (GPU) as a data processing accelerator and NVIDIA CUDA compiler 10.1 with the cuFFT library to speed-up the data analysis as described in our previous work [43]. The workstation was equipped with an Intel Core i7-7700K (4.2 GHz) CPU and 32 GB RAM. For data processing, we used NVIDIA GeForce GTX1080 with 8 GB of device memory.
We have measured the computation times of the implemented algorithms: PR-Doppler OCT, STdOCT, 2D 4FT FPU, and 3D 4FT FPU. When unwrapping the 3D data set with 2D 4FT FPU, the algorithm must be used in every B-scan separately. 3D 4FT FPU is called only once to unwrap the phase in the entire 3D data set. When a data processing accelerator is used (as in our system with GPU), there is a need to distinguish between the two processing platforms, i.e., the host and the device.The host is the CPU with RAM and the device consists of the GPU with its on-board memory. The bottleneck of the applications using GPU is the data transfer between the host memory and the device memory. Therefore, we present the computation time on GPU and the total processing time, including the raw OCT data transfer to the GPU memory, the GPU computations, and transfer of the results back to the host memory. Processing with and without data transfer was executed 10 times. The average computation times and standard deviations for protocols (a-c) from Table1 are listed in Table 2.

Comparison of PR-Doppler OCT and STdOCT
We compared the computation times and the axial flow velocity images obtained with PR-Doppler OCT and STdOCT to select the most appropriate method for imaging in the human eye. Figure 1 shows an example B-scan from a 3D data acquired in the human retina using protocol (c) from Table 1 with different combinations of Doppler OCT and 4FT FPU processing techniques. 3D 4FT FPU in combination with STdOCT rendered Doppler OCT images with visibly smaller noise as compared to other combinations and not affected by phase unwrapping artefacts. This observation is in agreement with our previous studies, where we showed that STdOCT offers higher SNR than PR-Doppler OCT [13,44] and, in another study, that 3D 4FT FPU is less likely to generate artifacts than 2D 4FT FPU [29]. It is interesting to relate these observations to the total computation times of the implemented algorithms presented in Table 2. Computations required for phase unwrapping of the entire 3D data sets with both phase unwrapping algorithms take less time than acquisition of the data. When combined with the PR-Doppler OCT, phase-unwrapped Doppler OCT images can be generated in times comparable with the data acquisition. However, computations required in the STdOCT method take half a minute when combined with the 3D 4FT FPU. This suggests possible uses of the two Doppler OCT methods. Real-time images are required during the OCT device alignment procedures, usually performed prior to the acquisition of the data to select the area to be imaged and to optimize the device settings. Typically, a series of 2D data is acquired for this "preview mode" and the images should be of acceptable quality if the highest quality images cannot be generated in real-time. A combination of PR-Doppler OCT and 2D 4FT FPU methods is well suited for this purpose. However, the 3D data acquired to aid the diagnosis of diseases should be processed with methods ensuring the highest signal-to-noise ratio and artifacts-free images. Combination of STdOCT and 3D 4FT FPU meets this requirement while still providing images 30 s after acquiring the 3D data set. Therefore, hereafter, all the results obtained in human eyes are presented for the combination of STdOCT and 3D 4FT FPU.

Flow phantom
The microfluidic device was positioned under the OCT microscope at an angle of ∼45°, as schematically illustrated in Fig. 2(a), to introduce a non-zero axial component of the flow velocity and to ensure visibility of the entire channel lumens in the C-scans. Experimental protocol (a) from Table 1 was used. The flow rate of the solution pumped through the microfluidic device was set to three different values, as listed in Table 3, at which a single phase wrap occurred in the Doppler OCT images, as shown in the first column of Fig. 2(b). The solution was injected into the large 300 × 100 µm channel and the flow then splits into three 100 × 100 µm channels. The green rectangles mark sets of N = 6 C-scans before and after trifurcation, which were used to obtain data presented in Table 3. The Roman numerals indicate the channels in the trifurcation. b) STdOCT Doppler C-scans selected from 3D data sets: one pre-and the second post-trifurcation. The two examples show data acquired for different volumetric flow rates, as labelled in the images. First column: wrapped phase. Second column: phase unwrapped with 3D 4FT FPU. The arrows indicate channels (I, II, III, IV), for which the flow rates are listed in Table 3 The flow rates calculated from the phase-unwrapped data in the regions of interest, marked by orange rectangles in Fig. 2(b), are listed in Table 3. For each flow rate set by the syringe pump (last column of Table 3), six adjacent C-scans were extracted from the data, in which the entire lumens of the channels were visible. The mean flow values were calculated for each channel. The uncertainty of the mean was estimated from t-Student distribution at 0.05 significance level. In a few cases phase unwrapping errors have occurred introducing errors to the calculated values (marked with crosses in Table 3). In the data with no phase unwrapping errors, the flow rates calculated for the large channel (channel IV) are equal to the flow rates set by the pump. This confirms that there are no offsets in the values calculated from the phase-unwrapped Doppler OCT data. As the net flow rate should remain constant, the sums of the flow rates in the trifurcated channels were also calculated. The flow rates calculated for the small channels (I, II, III) sum up to the flow rate in the large channel, except of those with identified phase unwrapping errors. This demonstrates that the phase-unwrapped data allows to reproduce the additive property of the volumetric flow rate.

Flow rates in vessels in the optic disc of the human eye
Doppler OCT imaging of five human volunteers (ID 1 to ID 5) was performed in the optic disc area. Six data sets according to protocol (b) from Table 1 were collected for each subject. Additionally, for subjects ID 1 to ID 3 six data sets according to protocol (c) from Table 1 were collected. Owing to the high blood flow rates and the steep angles between the vessels and the beam of light in this area of the eye, the axial component of the flow velocity vector often exceeds the velocity limit given by Eq. (4), and phase wrapping occurs. The flow rates were calculated in the STdOCT data unwrapped with 3D 4FT FPU method, in the manually fitted elliptical region of interest (ROI) for which we have calculated the flow rate according to Eq. (2). The minor axes of the ellipses were used as estimators of the vessel lumens. The mean values of the flow rates were computed from the six adjacent C-scans. The uncertainty of the mean estimation was calculated from the t-Student distribution at the significance level 0.05.

Vascular bifurcations
We have analyzed blood flow rates in five vascular bifurcations acquired from subjects ID 1 to ID 5 using protocol (b) from Table 1. In each subject we have selected a vascular bifurcation with vessels orientation in the 3D data sets as schematically illustrated in Fig. 3(a). Vessel III splits into two vessels I and II. Doppler OCT C-scans were selected at two different depths as listed in Table 4. Example C-scans proximal and distal to bifurcation are presented in Fig. 3(b).
In each of the five subjects six adjacent C-scans pre-and post-bifurcations were selected from which mean values were calculated and their uncertainties estimated from t-Student distribution at 0.05 significance level. Due to the scanning and the tilt of the vessels, branches I, II and III were imaged at different times. To minimize the change of the volumetric flow rate caused by the pulsatile blood flow, the C-scans were selected in such a way that the temporal separation of the sections of vessels I, II and III was less than 20 ms. Phase wrapping occurred in at least one of the branches. After phase-unwrapping, the volumetric flow rates were calculated for each of the branches. The results obtained from all five subjects are presented in Table 4. The additive property is fulfilled within the standard errors calculated from t-Student distribution at 0.05 significance level. Having vessels diameters d and from the elliptical fit we tested whether they fulfill the Murray's law [45], that links vessel diameter proximal to the bifurcation (vessel III) with the distal vessel diameters (vessels I and II) in the following way: d α In the original article α = 3 was found to be optimal value for vascular system that minimizes the power required to maintain and circulate the blood. In another derivation, where mass of the vascular system is minimized, α = 3 was found for laminar flows and α = 7/ 3 = 2.33(3) for fully turbulent flows [46]. In both cases the models assumed blood viscosity to be independent on vessel diameter, what is known to be inaccurate. The viscosity of blood depends not only on parameters such as hematocrit, pH, protein content, maturity of the cells (immature are more adhesive) but also depends on vessel diameter. It has the smallest value in the capillaries of diameter around 6-10 µm and increases in larger vessels [47]. In a recent study with almost 500 bifurcations measured with adaptive optics scanning laser ophthalmoscope values of α were observed to be higher in arteries than in veins and to be lower in smaller vessels than in larger vessels [48]. The authors reported α = 2.10 for arteries with diameters smaller than 50 µm and α = 2.74 for arteries diameters larger than 50 µm. For veins with lumen diameters smaller than 50 µm α = 1.69 and α = 2.11 for veins larger than 50 µm. The values calculated for bifurcations reported in the Table 4 for ID 1 (α = 2.17), ID 3 (α = 2.65) and ID 4 (α = 2.43) fall in these ranges, while values for ID 3 (α = 7.86) and ID 4 (α = 5.39) are outside. The latter is probably caused by the manual method of diameter estimations that is prone to errors if the visible vessel cross-section is not taken in the tubular part of vessel. The reason why the flow rate estimation is more resistant to this kind of error than the estimation of α is that the integration over erroneous ROI including parts of static tissue has no impact on the calculated flow rate as the static tissue component has zero mean velocity and does not contribute to the flow rate calculation.

Correlation between blood flow rates and vascular lumen diameters
Next, we have analyzed the relation between the blood flow rate and the diameters of vessel lumens in the images of subjects ID 1 to ID 3 using six data sets collected with protocol (c) from Table 1. We have identified a total of 600 vessels for which full lumen cross sections were visible in the Doppler OCT C-scans and estimated their diameters and flow rates using manually fitted ellipses as described before and shown in Fig. 4. The diameters were in the range of 12-150 µm. We have excluded capillary vessels with diameters smaller than 12 µm, which is the lateral imaging resolution limit of the OCT system. We have classified the identified vessels as arterioles and veins using Doppler OCT information about flow velocity direction in the center of the optic disc. Next the branching of the vessels was followed in the OCT angiography images as well as in the structural OCT images shown in Fig. 5, generated from the same data sets as Doppler OCT images, where we were able to identify which vessel originated from the main arteriole and which is a vein. Vessels with lumen diameters smaller than 20 µm were classified as capillaries.
In Fig. 6 we have plotted the values of the flow rates vs. lumen diameters in the phase-wrapped and phase-unwrapped data in a linear scale graph for subject ID 1. The blood flow rates calculated in the vessels with phase-wrapping, presented as red dots in Fig. 6(a), seem to be randomly distributed with values below ∼11 µl/min suggesting no relation between the vessels diameters and blood flow rates. In the vessels where no phase-wrapping occurred, presented as blue crosses in Fig. 6(a), the data shows a power dependence Q ∼ d α , known from literature [36], between the blood flow rate Q and the lumen diameter d. After performing the phase-unwrapping with the 3D 4FT FPU all the detected vessels show that dependence between the vascular lumens and the blood flow rates. The fit with a power curve was performed using a weighted least-squares method and is presented in Fig. 6(b). The correlation parameter was found to be α = 1.8 ± 0.2 for veins and α = 2.0 ± 0.2 for arterioles. The results obtained from the 600 vessels identified in 18 data sets of three subjects are presented in Fig. 7 as log-log plots of the flow rates versus the lumen diameters. The calculated flow rates were in the range of 0.16-25.14 µl/min for the veins and 0.15-18.32 µl/min for the arterioles. The plots indicate a positive linear correlation in the log-log scale. The slopes of the linear fit (or correlation coefficients α in the power curve fits) were found to be equal to α = 1.8 ± 0.2, α = 1.9 ± 0.2, α = 2.0 ± 0.2 in case of veins and α = 2.0 ± 0.2 α = 1.8 ± 0.2, α = 2.2 ± 0.2 in case of arterioles for subjects ID 1 to ID 3, respectively. In case when the capillaries with diameters less than 20 µm are not taken into account in the fitting procedure, then the above values increase to α = 2.0 ± 0.2, α = 1.9 ± 0.2, α = 2.0 ± 0.2 in case of veins and α = 2.4 ± 0.2, α = 1.9 ± 0.2, α = 2.2 ± 0.2 in case of arterioles. This may suggest that either systematic errors appear in vessel diameter estimation and/or the flow mechanics is different in vessels with various diameters as mentioned earlier [47,48]. The above values of slopes are in agreement with values found by other researchers. For example Wang et al. [31] reported slopes in the range 1.52 − 2.54 with mean slope 1.97 ± 0.40. Higher slopes, closer to Murray's value of three [45] can also be found in the literature. Riva   a Circular standard deviation of the phase distribution defines the noise level of input data [29]. et al. [48] in their study on retinal bifurcations. If the assumptions made in the derivation of Murray's law are correct [45], then the slope of the linear fit and the power coefficient in the law governing the diameters of bifurcating vessels have the same value. Therefore, as mentioned before in Section 3.3.1 only large arteries with diameters above 50 µm show slopes close to three (2.74), while all the other vessels have slopes around two. The results presented in Fig. 7 are in agreement with the results presented by other studies and show that phase unwrapping is a necessary step in the flow estimation. The 3D 4FT FPU approach to phase unwrapping is a good candidate to become a standard step in the Doppler OCT data processing pipeline due to its robustness and high computational efficiency. It has to be noted, though, that the weakest point in the current study is the procedure of vessel diameter estimation. First, the precision of this estimation was limited by the transverse imaging resolution of the OCT device equal to 12 µm. Second, the estimation of the lumen diameters was done in C-scans by manually fitting ellipses and using their minor axes as lumen diameters. In this procedure it is not always obvious how the ROI should be placed, especially in the proximities of bifurcations, where the cross-sections of the vessels are not always ellipsoidal, what can be observed in Fig. 4. Although this manual method is prone to errors the general conclusion is visible -the retinal blood flow in the vessels obeys the power dependence on vessel lumen diameter. Development and optimization of procedures yielding more precise values of vessel lumen diameter is outside of the scope of this paper but will be addressed in our future works.  (Table 4.). Top: selected cross-sectional images of the wrapped phase. Bottom: outcomes of the 3D 4FT FPU algorithm. The arrows with Roman numerals indicate the selected vessels. The subscript "W" refers to the wrapped phase, while the subscript "U" refers to the same vessels after unwrapping. The figure insets show magnified regions of the indicated vessels. I U , III U : errors of the phase unwrapping (phase has remained wrapped). Arrows II U , IV U : phase-wrap-free vessel (diameter: 18 µm, 31 µm). Arrow V U : error-free phase unwrapping of the flow for a vessel with diameter of 79 µm. The insets for vessel II U present five zoomed regions of adjacent C-scans. The indicated vessels are plotted in Fig. 7    Log-log plots of the blood flow rate versus the vessel lumen diameters. First row: venous blood flow. Second row: arterial blood flow. The green marks indicate capillaries with a diameter smaller than 20 µm. Each column shows combined results from 6 data sets obtained in one of the three subjects (ID 1 to ID 3, Table 4). Solid lines: linear fit (the equations obtained by the least-squares fit are given on the top of each plot). Dashed lines: 95% confidence interval of the fit. Example vessels I U , II U , IV U , V U , and VI U presented in Fig. 4 are indicated with magenta outlines.

Conclusions
Phase unwrapping is a procedure which should be implemented in the Doppler OCT algorithms for calculation of the volumetric blood flow rates in the wide range of flow velocities encountered in the vessels of the human retina and optic nerve head. We have demonstrated that the 3D 4FT FPU algorithm applied to Doppler OCT data gives reliable measurement of the volumetric flow rates. Since there are no "gold standards" for the volumetric blood flow rates in the vessels of the eye fundus to which we could relate our results, we have performed experiments in a flow phantom and tested relative relations between blood flow rates in the eyes of volunteers. The conclusion from the flow phantom experiments is that not only the measured flow rates fulfill the additive property of split-flows but also, they are in an agreement with the flow rates set in the microfluidic pump. The analysis of the blood flow rates in the vascular bifurcations of human subjects demonstrated that the relations between flow rates, measured as the additive property of flow, are preserved. Since there are no offsets between induced and measured flow rates in the flow phantom experiments, it is unlikely that such offsets should be introduced by our methods to the in vivo Doppler OCT imaging of the eye. Therefore, the results of the volumetric blood flow rate measurements in the 600 vessels identified in 18 data sets acquired in the eyes of three volunteers may be considered as reasonable estimation of the actual blood flow rates.
As a demonstration of the necessity of the phase-unwrapping in the Doppler OCT, we have analyzed the relation between the blood flow rates and the vessel lumen diameters. We have observed that the phase-wrapping causes large dispersion of flow rates especially for vessels with large diameters. The flow rates included unexpectedly low values for large vessels. Only data free from phase-wrapping (either inherently or after application of the phase-unwrapping methods) reveals the positive power dependence between the blood flow rate and vascular lumen diameters in accordance with theoretical predictions and reports of other researchers.
We also introduced the 4FT FPU methods into the OCT processing pipeline implemented in GPU showing that phase-unwrapped phase-resolved Doppler OCT can be processed in real-time, and phase-unwrapped joint spectral and time domain OCT images can be generated in less than a minute (Table 2). This is an important result which shows that phase-unwrapped Doppler OCT methods can be used in clinical practice.