Blood flow velocity vector field reconstruction from dual-beam bidirectional Doppler OCT measurements in retinal veins.

In this paper, we demonstrate the possibility to reconstruct the actual blood flow velocity vector field in retinal microvessels from dual-beam bidirectional Doppler optical coherence tomography measurements. First, for a better understanding of measured phase patterns, several flow situations were simulated on the basis of the known dual beam measurement geometry. We were able to extract the vector field parameters that determine the measured phase pattern, allowing for the development of an algorithm to reconstruct the velocity vector field from measured phase data. In a next step, measurements were performed at a straight vessel section and at a venous convergence; the obtained phase data were evaluated by means of the new approach. For the straight vessel section, the reconstructed flow velocity vector field yielded a parabolic flow. For the venous convergence, however, the reconstructed vector field deviated from a parabolic profile, but was in very good accordance with the simulated vector field for the given vessel geometry. The proposed algorithm allows predictions of the velocity vector field. Moreover, the algorithm is also sensitive to directional changes of the flow velocity as small as <1°, thereby offering insight in the flow characteristics of the non-Newtonian fluid blood in microvessels.


Introduction
Optical coherence tomography (OCT) is a technique that uses low coherence interferometry to obtain tomograms with a resolution in the micrometer range. In ophthalmology, OCT has revolutionized imaging and has been widely used for diagnosis and therapy monitoring [1,2]. In recent years, OCT has also been extended towards functional applications such as Doppler OCT [3], which combines the principles of laser Doppler techniques [4] with those of OCT imaging. In Fourier-Domain Doppler OCT, information on the blood velocity is obtained by phase sensitive detection of the complex OCT signal i.e. calculating the phase shift by subtracting adjacent A-scans [5,6], in the following referred to as phase. To measure absolute velocities in vivo, however, the problem of the a priori unknown Doppler angle has to be solved. One solution is to employ a dual-beam bidirectional Doppler system; this allows for absolute blood velocity measurements independent of the Doppler angle [7][8][9].
The retinal microcirculation is the only vascular bed that offers a direct view on micro vessels. A wide variety of studies have focused on retinal vascular calibers because they are easily accessible. Small retinal arterial diameters have been shown to be associated with cardiovascular disease [ In the present paper, we propose a technique to extract the velocity profiles at retinal bifurcations based on dual-beam bidirectional Doppler OCT. To the best of our knowledge, this is the first approach to measure 3-dimensional blood flow velocity profiles at bifurcations of microvessels in the order of 100-200 μm in any vascular bed of the human body. The technique may have considerable potential for studying patients with systemic, cerebral and/or ocular vascular disease. Furthermore, studying velocity profiles should improve the understanding of the flow characteristics of the non-Newtonian fluid blood in microvessels. This may significantly increase our understanding of ocular perfusion in diseased states.

Dual beam Doppler Fourier-domain OCT
The measurements were performed with a dual-beam bidirectional Doppler OCT system that is depicted in Fig. 1 and described in detail elsewhere [9]. The system uses a superluminescent diode (SLD, Superlum, Carrigtwohill, Cork, Ireland) (λ = 838 nm, Δλ = 54 nm) as light source. The light of the SLD is split into two probe beams separated by their polarization properties and impinging onto the vessels of interest under a known angle. The spectral signal of these two channels is recorded by two identical spectrometers; thus, the phase shifts caused by the moving red blood cells (RBCs) in the sample is recorded from two different angles. To obtain a live image of the fundus, the OCT system is coupled to a modified fundus camera (FF450plus; Carl Zeiss Meditec AG, Jena, Germany) by means of a dichroic mirror. This fundus camera and a charge coupled device (CCD) camera with hardand software units for vessel analysis (Dynamic Vessel Analyzer, DVA, Imedos, Jena, Germany) are used for the examination of the vessel diameters at the measurement position. Moreover, the fundus view allows the operators to choose the desired measurement region at the ocular fundus and enables the accurate positioning and focusing of the probe beams.

Measurement technique
The measurement protocol was approved by the Ethics Committee of the Medical University of Vienna. OCT and vessel diameter measurements shown in this paper were performed on two healthy subjects, both aged 28. The nature of the measurements was explained to both subjects and they gave written consent to participate. The protocol of the study complied with the standards of the Declaration of Helsinki.
To determine the flow velocity field at a venous bifurcation, scans were performed at the convergence of interest, such as the exemplary vascular section shown in Fig. 2(a). After choosing a suitable position, great care was taken to focus both beams onto exactly this position in the fundus. This was done in an effort to guarantee a complete knowledge of the measurement geometry, i.e. the separation angle between the two probe beams, necessary for reconstruction of the blood flow velocity field from dual beam Doppler FD-OCT measurements.
To allow for a fast recording, the number of A-scans per B-scan was limited to 500. Depending on the expected blood velocity at the measuring site, the CCD cameras' integration time was set to a value between 25 µs and 60 µs per A-scan. To ensure adequate oversampling, a relatively small scanning width was chosen, usually about 6 vessel diameters (approx. 1000 µm). It would be unfavorable to choose a smaller scanning width, as the latter had to be large enough to ensure that the vessel stays in the scanning area over the whole recording time despite the inevitable subject motion. For the results shown in this paper, between 60 and 300 frames were recorded for both channels at each selected measuring position in the fundus.

Preparation of the phase images
To obtain the phase images, all frames were bulk motion corrected [17,18] and automatically registered to a user selected reference frame. For each frame and the reference frame, a metric was calculated to identify potentially falsely aligned or motion corrupted frames. This allowed the rejection of frames for which the metric exceeded an empirically defined threshold. The portion of these frames that deviated from the reference frame was in the range from 0 to 5 percent of the whole data set. After these first post-processing steps, all remaining phase frames were assumed to originate from the same position within the fundus. In a next step, the phase frames were averaged separately for channel 0 and channel 1. To lower the amount of data to be handled and allow a comparison of the two averaged phase frames, the frames were cropped to a rectangular area containing only the vessels. The two resulting matrices represent the two averaged phase frames; an example of these matrices is shown in Fig. 3.

Simulation of phase profiles
As can be seen in Fig. 3, the phase profiles of the two channels do not equal although they were measured at the same fundus position. Another salient feature that merits attention is the sign change of the phase within the vessel area. To understand the nature of the averaged phase profiles, the scanning geometry of the system has to be considered.
The motion of scattering particles, in this case the RBCs in the vessel, causes a phase shift in the back reflected light. However, the measured phase values represent only the component of the velocity vector pointing in direction of the incident probe beam. Figure 4(a) shows the 3D scanning geometry for one of the two channels: The yellow line illustrates one of the two probe beams, which are separated by an angle α and both of which lie in the y-z-plane. The bold red arrow represents the velocity vector v  of the erythrocytes flowing within the vessel, which is tilted by an angle γ with respect to the plane orthogonal to the optical axis of the system (blue) and by an angle β with respect to the y-z-plane, i.e. the plane spanned by the probe beams. To obtain the component of the velocity vector that actually causes the measured phase pattern, one must first calculate the orthogonal projection of the velocity vector onto the plane spanned by the scanning beams (the y-z-plane in Fig. 4(a)). This yields the dashed red arrow v′  .
Hence, || v′ can be calculated using the following formula: As a consequence, the measured phase profiles cannot be converted into velocity flow profiles if the angles β and γ are unknown. Conventionally, when measuring the mean flow velocity within a vessel cross-section with dual-beam Doppler FD-OCT, these angle dependencies are overcome by measuring at vessel sites where the vessel forms a straight line and lies flat in the fundus. This allows two assumptions that facilitate calculations: firstly that the angle β (angle between the plane spanned by the two scanning beams and the blood flow velocity vector) remains constant over the vessel cross section and can simply be deduced from the fundus image or an en-face image taken while measuring, and secondly that the angle 0 For measurements at venous convergences as described above, these two assumptions are no longer valid. A constant angle β over the whole vessel width is unlikely at the convergence and γ must change since the two vessels before and the one after the convergence have different diameters. The influence of changes in flow direction (i.e. changing angles β and γ ) on the measured phase becomes clear upon closer examination of Eq. (3): if scanning is performed perpendicularly to a vessel, β is close to 0° and its influence on the flow as the argument of the cosine is low. However, the argument of the second cosine -( 90 2  As expected, a parabolic velocity vector length distribution causes a parabolic phase profile when γ is constant over the vessel depth. At convergences, however, the diameter of the vessel after the convergence is larger than the diameter of each of the two joining vessels, which leads to a divergent flow direction and accordingly to a non-constant angle γ . To simulate the phase profile at a junction where the vessel geometry widens, the angle γ of the velocity vectors was adjusted accordingly: the vector at the vessel bottom was set to −10°. The tilt γ of the remaining velocity vectors over the vessel depth was increased linearly, reaching 10° for the vector at the top of the vessel and giving the situation depicted in Fig.  6(a). Obviously, an increase in vessel height, which leads to diverging velocity vectors, causes a sign change in the phase profile (see Fig. 6(b)). This corresponds exactly to the observation that was made in the in vivo measurements at the venous convergence (see Fig.  3). The observed phase shift is also dependent on the angle between the probe beam and the vessel. Changing the angle under which the scanning beam impinges onto the vessel from 2.75° (as in Fig. 6) to −2.75° (as shown in Fig. 7), results in the situation in the second channel of the dual-beam OCT system (see Fig. 3(b)): evidently, the zero crossing of the measured phase now takes place at a different depth position within the vessel. This is due to the fact that the phase is zero when the scanning beam impinges perpendicularly to the velocity vector. The incident angle of the scanning beams is different for each channel; consequently, a velocity vector with a distinct angle γ at a different distinct vessel depth is perpendicular to one or the other scanning beam and causes a phase of zero.

Analysis of the measured phase images and estimation of the blood flow velocity field
The previous section illustrated the influence of changes in γ on phase profiles in systems zdirection measured using two scanning beams separated by an angle α; thus the unexpected phase changes in the measured data could be explained. With this information, it was possible to reverse the process, i.e. to estimate the blood flow velocity fields that caused specific measured phase profiles. This is necessary because knowing phase profiles only reveals little information on the actual blood flow velocity field at the measurement site.
To gain insight into the velocity field, we took the following approach (see Fig. 8): in the images from both channels, corresponding areas of adjacent A-scans (usually of a width of about 5 pixels and ranging from the top to the bottom of the matrices) were chosen (marked by black bars in Fig. 8(a)) and averaged in x-direction for smoothing. This yielded onedimensional arrays of phase values for each channel for a selected lateral range (Fig. 8(b)). We then searched for a velocity field that could cause the measured phase depth profiles. To do so, a MATLAB® script (MATLAB® R2013b, The MathWorks Inc., Natick, Massachusetts) was written which simulated a velocity vector field, calculated the phase profile which this velocity vector field would cause in channel 0 and channel 1 and compared these simulated phase values with the measured phase data set. To find the most suitable velocity vector field, the program varied different coefficients as outlined in detail below.
First, a list of scalars i r was generated. Each of the values i r corresponds to the length of a velocity vector at equidistantly spaced depth positions i k (whereby the vessel is defined to start at 10 i k = − and end at 10 i k = ) and was calculated with the parameters b and τ as follows: The initial value for b was set to 2, which means that the vector lengths are distributed parabolically. Due to the flow characteristics of blood and the merging situation at convergences, a perfectly parabolic velocity flow profile is unlikely and so variation of b by the algorithm is allowed.
The coefficient τ is a scaling factor for the values of the scalars (which, in turn, correspond to the lengths of the velocity vectors). Since the length of a vector influences the magnitude of the measured phase, a variation of τ by the algorithm is intended.
Thereafter, the γ -distribution over the vessel depth had to be simulated. To accomplish this, the algorithm could change the slope s γ of the γ -distribution and add an offset off γ : { }  Figure 9(a) shows the best fits (green) for the averaged phase depth profiles (blue) shown in Fig. 8(b). Figure 9(b) shows the velocity depth profile the method yields for the area between the black bars in Fig. 8  To gain an insight into the blood flow vector field distribution over the whole vessel area, the algorithm was applied to equidistant lateral (x-) positions (x start = 20, x end = 80) over the vessel width. The determined coefficients are plotted in Fig. 10.   The method was also applied to the measurements performed at the venous convergence (Fig. 12). In Fig. 13, the results of the calculation for the averaged depth scan shown in Fig. 12 are displayed. The result of the variation of the coefficients is presented in Fig. 13(a). The blue line depicts the measured phase while the green line shows the phase of the simulated vector field which has the smallest divergence from the measured phase pattern. Clearly, whereas the overall flow area decreased as compared to the situation before the junction site (vessel diameters before junction: 134.6 ± 3.3 and 126.7 ± 3.5 µm; diameter after junction: 163.2 ± 3.1 µm), the axial extent of this area (z-direction in Fig. 13(b)) increases. This can also be substantiated by numbers: as the phase values chosen for the evaluation were taken from the left side of the phase image (black bars in Fig. 12(a)), the flow at this position after the convergence originated from the left of the two vessels. Before the convergence, this vessel had a diameter of 126.7 µm while the vessel diameter after the convergence was measured to be 163.2 µm. This means that the flow must diverge in z-direction with a strong gradient. Figure 14 shows a supposable vessel wall geometry (thick red lines) overlaid over the calculated velocity vector field. Under the assumption of a linearly increasing vessel diameter from 127 ≈ µm to 163 ≈ µm, the vessel segment over which this increase takes place is 112 ≈ µm. To obtain the velocity vector field over the whole vessel area, the algorithm was applied at equidistant segments in x-direction over the vessel width (step size 20 pixels). The determined coefficients are plotted in Fig. 15.  Figure 16(b) shows an overlay of the calculated velocity vector field at a distinct xposition with the measured phase profile of one channel. In Fig. 16(a), the phase difference of the two channels is plotted.

Discussion
The introduction of OCT has revolutionized ophthalmology; in the meantime, OCT has gained an important role in the diagnosis and treatment monitoring of ocular pathologies. Doppler OCT is an attractive extension of the technique [3] allowing for the quantification of blood flow [9, 11, 14, 21] as well as for non-invasive angiography [22][23][24][25]. The present work that is providing information about the velocity profiles at retinal bifurcations may add to this wide array of OCT applications, since it gives insight in the actual flow conditions at the measuring site and reveals details of the flow characteristics of blood in microvessels. This is of interest because a variety of common eye diseases such as glaucoma [26, 27], age-related macular degeneration [28] and diabetes [29] is associated with ocular perfusion abnormalities. In addition, the retinal microvasculature is an important biomarker to gain insight into systemic cardiovascular diseases [30,31].
In this paper, we report on a new algorithm to extract velocity profiles in microvessels which was applied to two different flow situations that were measured by means of dual-beam bidirectional Doppler OCT. Due to the geometrical complexity it is not feasible to reconstruct the blood flow velocity vector field from one of the channels presented in Fig. 3. The use of two linear independent probing beams is necessary to obtain a unique flow vector field. First, a straight vein (diam. 98 µm) was measured for which the algorithm yields the coefficients shown in Fig. 10: off γ is almost constantly 4° and s γ is almost zero over the whole vessel diameter. The slight deviations from zero at the vessel borders might stem from a lower signal to noise ratio (with respect to the phase signal) near the vessel wall where the phase values are small. The coefficient τ, which is proportional to the mean blood flow velocity at a position x in the vessel, shows an almost parabolic, slightly flattened velocity vector distribution. The coefficient b, which lies between 2.2 and 2.9, indicates a slightly flattened velocity profile as well.
The second examined flow situation was a venous convergence. The measured phase profiles looked different and showed sign changes within the phase pattern over the vessel area. With the applied algorithm, it was possible to understand these sign changes and to estimate the velocity field underlying the measured phase profiles. The resulting coefficients are shown in Fig. 15. The offset angle off γ changes from about −2° to −5° over the convergence width. This indicates that the two source vessels join the main vessel under different angles with respect to the optical axis of the system. The coefficient s γ suggests an increasing diameter of the flow pattern in z-direction for the left vessel and equally, but less pronounced, an increase of diameter in z-direction for the right vessel. These results are in good agreement with the measured vessel diameters before and after the convergence. The coefficient τ shows a slightly flattened velocity profile over the vessel width and coefficient b also indicates flattened profiles in z-direction. This seems plausible when considering the merging flows and the changing flow geometry in z-direction. To show the advantage of a detailed analysis of the measured phase profiles, an overlay of the measured phase and the calculated velocity vector field for the convergence is presented in Fig. 16(b). The divergence of the velocity vectors over the vessel depth can clearly be seen. When comparing the outcomes of the algorithm (Fig. 16(b)) to the usual approach of subtracting one phase image from the other (Fig. 16(a)), it is obvious that the correct vessel geometry as well as the full direction of the blood flow vectors can only be gained by applying the reconstruction algorithm.
In the present model, we did not take into account changes in angle β between the plane spanned by the two scanning beams and the blood flow velocity vector. In the case of the straight vessel segment, this angle can easily be determined from the fundus image. In the case of the convergence, this question is more complicated. However, its influence on the measured phase as the argument of the cosine is low. Since the calculations are performed on averaged depth scans around the same lateral (x-) position, it is unlikely that within this relatively narrow interval the angle β changes significantly over the vessel depth. As a consequence, the influence of changes in β over averaged depth scans on the calculations is considered to be low.
In principle, vessels that are measured at the ocular fundus by means of Doppler OCT are arterioles and venules with a diameter between 50 and 200 micrometers. The measurements presented in the manuscript in hand were performed in retinal veins. Basically, the algorithm is also applicable to data that stem from retinal arterial sections. However, due to the high pulsatility in arteries, a pulse triggered OCT recording [32] is required to conduct the frame averaging that is necessary for the reconstruction of the velocity vector field.
Indeed, only very few velocity profile data are available for convergent venous vessels of the same size as those presented here from the human retina. Recently, Leble and associates presented en face in vitro data using a microchannel containing both a diverging and a converging bifurcation. In fact, the authors also obtained non-parabolic velocity profiles when the two vessels were close to the point of convergence [33]. Trasischker and associates introduced an approach using three probe beams allowing to reconstruct spatial orientation and magnitude of the blood flow in a vessel without prior knowledge of the measuring geometry [34].
To the best of our knowledge, the data presented in this manuscript are the first on blood velocity profiles in human vessels offering a reconstruction of the full three-dimensional velocity vector field. Previous data have only shown one-dimensional profiles using scanning laser Doppler velocimetry [35] or microparticle tracking velocimetry [33], where the first lacks depth information and the second does not provide information about movements parallel to the probe beam within the sampled vessel. Previously published measurements with standard Doppler OCT showed relative profiles that have not been corrected for the axial angle [5].

Conclusion and outlook
This paper presents a new analysis method for dual-beam bidirectional Doppler OCT measurements. Two flow situations (straight vessel, convergent vessels) were firstly simulated, and, in a next step, measured and analyzed with a newly developed mathematical algorithm. The approach takes the full geometrical complexity when assessing the velocity information recorded by means of dual-beam Doppler OCT into account. The resulting data is in good agreement with the expected values and allows a reconstruction of the actual blood flow velocity vector field in retinal vessels, which, to the best of our knowledge, has not been shown before. The results obtained in healthy subjects seem promising and the analysis method offers a high potential for examining retinal blood flow abnormalities in patients with ocular diseases and gaining an insight into the flow characteristics of the non-Newtonian fluid blood at vessel diameters down to 100 µm or even smaller.