2.1.

Theory of En Face Method for TRBF Measurement

The theory of flow calculation using Doppler OCT and en face plane has been previously described.^34,36 Briefly (Fig. 1), the Doppler phase shift is proportional to the axial component of velocity, or $D\propto V\cos (\alpha )$, where $D$ is Doppler phase shift, $V$ is velocity, and $\alpha$ is the angle between the axial component and vessel. The true vessel area, which is perpendicular to the velocity vector in vessels, is equal to the product of vessel area in the en face plane and $\cos (\alpha )$, or $S=S_e\cos (\alpha )$. Here, $S$ is the area of a vessel in the direction perpendicular to the vessel segment and $S_e$ is the area of the vessel in the en face plane. When we do the integration of V in $S$, the two $\cos (\alpha )$ components cancel each other. Thus, blood flow in a vessel is proportional to the integration of Doppler phase shift in the vessel region and en face plane.

Fig. 1

Illustration of en face Doppler technique for blood flow calculation in a vessel.

2.2.

Doppler OCT and Scan Pattern

The OCT system that we used in this study was RTVue-XR spectral OCT. RTVue-XR has a 70-kHz scan rate, a laser source with central wavelength of 840 nm and $5\text{-}\mu \mathrm{m}$ depth resolution. The phase wrapping velocity limit of RTVue-XR is $11.1\mathrm{mm}/\mathrm{s}$ using the formula $V_{\text{max}}=\lambda /(4n∆T)$ and wavelength $\lambda =840\mathrm{nm}$, refraction index $n=1.36$ and axial scan (A-scan) time interval $∆T=14.3\mu \mathrm{s}$.

The highest blood velocity in the central retinal veins has been quoted as up to $57\mathrm{mm}/\mathrm{s}$ using ultrasound.³⁸ Because the central retinal vein lies along the axial direction, the axial speed would exceed multiple phase wrapping limit of the 70-kHz OCT system. On the other hand, the average speed in peripapillary veins is 15 to $16.8\mathrm{mm}/\mathrm{s}$ as assessed by Doppler OCT.^6,39 In this area, the vessels are closer to perpendicular to the axial direction. Thus, the axial speed in peripapillary veins might be too small to be measured precisely. Therefore, we targeted vein branches near the disc boundary, in which speeds are between peripapillary veins and center retinal veins. The axial speed in these branches can be reliably measured when the angle is appropriate in some en face planes.

The scan pattern for TRBF measurement contained five volume scans covering $2\times 2\mathrm{mm}$ area around the optic disc. The volume scan contained 80 B-scans and each B-scan contained 500 axial scans. The scan depth was 2.3 mm. The volume scan was repeated five times in the scan pattern and took 3 s to obtain (Fig. 2).

Fig. 2

Scan pattern used for Doppler optical coherence tomography (OCT). The raster scan, or three-dimensional (3-D) volume scan, covered a $2\times 2\mathrm{mm}$ area around the central retinal vessels and was repeated five times in 3 s.

2.3.

Doppler Phase Shift Calculation and Bulk Motion Removal

The Doppler phase shift was calculated based on the phase-resolved technique.^18,23 With a faster scan rate, the noise is more than last generation of RTVue system.⁴⁰ In order to reduce the noise in OCT image, we applied the split spectrum method, which applied a filter bank to divide the full-spectrum fringe into different bands.^41,42 The key of split spectrum method is to average Doppler phase shift in all bands and reduce the noise. The cost is lower depth resolution due to the smaller spectral bandwidth in each band.

First, we generated the split spectrum by passing full-spectrum fringe through filter bank. Then, we applied the fast Fourier transform (FFT) transform to each band separately [Eq. (1)]:

Here, $S(x,k)$ is the spectrum obtained for $x$’th A-scan in k-space, $f(j,k)$ is the $j$’th band pass filter, and $s(x,j,z)$ is the complex signal of $j$’th band of the $x$’th A-scan.

Doppler phase shift was then averaged in all bands and neighboring A-scans [Eq. (2)]:

Here, $D(x,z)$ is the Doppler phase shift of the $x$’th A-scan at depth position $z$, ARG is a function to give the angle of a complex number, $N$ is the number of A-scans in the neighbor, $M$ is the number of bands, and $\mathrm{\Delta }T$ is the time difference between $x$’th and ($x+1$)’th A-scans. By averaging the $N$ A-scans and $M$ bands, the noise of the Doppler OCT image was reduced.

Bulk motion due to eye movement along the depth direction was calculated and removed using the averaged Doppler phase shift of the static retinal tissue as the bulk motion induced phase shift. Here, a revised histogram estimation method for bulk motion was used.⁴³ The original method calculated the histogram of the Doppler phase shift in static retinal tissue for an A-scan and used the value of the maximum bias bulk motion. In regions with less effective data points such as regions with large vessels and a thin retina, the method may not be precise. As each B-scan was acquired in less than 10 ms and the bulk motion changed gradually along the B-scan, a weighted polynomial fitting was applied on the bulk motion estimated by original method. To reduce error, weights were set to 0 for A-scans with less effective points.

2.4.

Doppler Phase Shift Calibration Experiment

Beyond the phase wrapping limit, fringe washout reduces the OCT signal strength and increases phase noise. Therefore, phase-unwrapping algorithms often do not recover the full Doppler shift signal and have a tendency to underestimate flow speed. It has been reported that the relationship between Doppler shift obtained from the phase-unwrapping method and the actual speed was not perfectly linear.⁴⁴ On the other hand, when the space interval between two consecutive axial scans increases, the correlation between two axial scans becomes weak and noise of the estimated Doppler phase shift increases. Therefore, the summation of Doppler phase shift in the vessel area will be smaller than the actual Doppler phase shift produced by blood flow. It was called the step-size factor and needs to be compensated.²³ In this study, we compensated for both bias by empirical calibration using a physical flow phantom.

A syringe was filled with fresh bovine blood (Animal Technologies, Inc., Tyler, Texas) and then pumped by a syringe pump (11 plus, Harvard Apparatus, Holliston, Massachusetts). The blood was pumped with a constant speed through a plastic tube with a $200\mu \mathrm{m}$ inner diameter [Fig. 3(a)]. The tube was fixed in front of the OCT scanner. The angle between the tube and laser beam could be adjusted using an articulating base [SL20, Thorlabs, Inc., Newton, New Jersey, Fig. 3(b)]. An achromatic doublet (AC254-030-B, Thorlabs, Inc.) was placed between the OCT scanner and tube in order to focus the laser beam on to the tube. The same scan pattern described earlier for TRBF measurement was used to obtain the 3-D image of the tube. In the experiment, the pump flow rate was set to be $60\mu \mathrm{l}/\min$, which correspond to an average velocity of $31.8\mathrm{mm}/\mathrm{s}$ in the $200\text{-}\mu \mathrm{m}$ tube. The speed was specifically chosen to simulate the blood flow rate in retina vessels inside optic disc. We first adjusted the vertical and horizontal position of the tube to center it with laser beam [Fig. 3(c)]. Then, we changed the incident angle from 0 to 18 deg with steps of 3 deg. The angle was limited to 18 deg because the reflectance signal is too weak due to fringe washout when the angle is larger than 18 deg.

Fig. 3

Flow phantom to calibrate Doppler phase shift and speed. (a) Flow phantom and OCT scanner; (b) articulating bases for rotation; (c) centration of OCT scan on tube of flow phantom (left: en face projection, right: a horizontal OCT scan at the center of 3-D scan pattern). (d) Parabolic distribution of blood flow in the tube. The velocity at the center of tube is about two times the overall average speed $V_{\mathrm{avg}}$; and (e) transverse position of three frames (or B-scans) and the corresponding cross-section image.

We used the same TRBF scan pattern for the phantom test in order to keep the same A-scan space interval as in vivo. Therefore, we need not consider compensation for the step-size factor after the calibration. The Doppler phase shift was calculated using the method described in the last section. Due the difference of reflectance index between air, tube, and blood, there is significant shape distortion at the edges of the tube where the incident angle is much larger than zero. To avoid this distortion, we only used the signal in the center of the tube, where the incident angle is close to 0 and the distortion is negligible. The flow in the tube followed a parabolic distribution [Fig. 3(d)] along the tube radial direction in theory. Therefore, the Doppler shift in the center of tube was larger than at the edges. We recalculated the theoretical average speed in the central region with radius $<20\mu \mathrm{m}$ to be $61.4\mathrm{mm}/\mathrm{s}$. The average velocity was obtained by averaging the Doppler shift in the center part of tube in the Doppler OCT image. Image segmentation was applied on reflectance image to find the top position of tube on each B-scan. Then, we estimated the beam incident angle based on the 3-D orientation of the tube [Fig. 3(e)]. The phase unwrapping algorithm is used when it is necessary. The axial speed was then calculated based on the Doppler shift and estimated incident angle.

2.5.

Pulsatility in Vessels

Pulsatility during the cardiac cycle makes the blood flow in the retina variable. To obtain repeatable TRBF, we needed multiple volumetric scans to match different time points in a cardiac cycle. It was difficult to obtain a high volumetric scan rate in the system used in this study. In order to cover the entire optic disc, the scan size needed to be about $2\times 2\mathrm{mm}$. We needed to limit the A-scan interval to ensure the correlation of neighboring A-scans and also to limit the B-scan interval to ensure that we included middle size vessels ($D>50\mu \mathrm{m}$). Based on these constraints, we chose a scan pattern with 80 B-scan $\times 500$ A-scans, covering a $2\times 2\mathrm{mm}$ area. The volumetric scan rate is 1.7 Hz.

The pulsatility of an artery cannot be compensated with such a volumetric scan rate. We compared the pulsatility profile of a retinal artery and vein using a rapidly repeating circular scan. In this experiment, eyes were scanned by repeating circular scans with 1.5-mm diameter at 30 Hz. In each cross-sectional image, the Doppler phase shift in the artery and vein was summated in the vessel area and normalized by its own mean. The start position of each cardiac cycle was detected using the maximum change in artery flow profile. Then, the profile of all cardiac cycles was aligned and averaged. By this method, we get an average pulsatility profile of artery and vein in one cardiac cycle (Fig. 4). We found that the variance in an artery was quite large compared to its mean ($\text{coefficient of variation}=0.26$), whereas the variance in a vein was relatively small compared to its mean ($\text{coefficient of variation}=0.05$). Thus, the pulsatility variation can be compensated for if we measure only venous flow and average enough number of volumes. The slower flow speed in veins also caused less phase wrapping than in arteries.

Fig. 4

Pulsatility profile of artery and vein.

2.6.

Vessel Detection and Phase Unwrapping

In order to detect veins, we developed an automated algorithm [Fig. 5(a)]. First, the retina was detected based on the reflectance image [Fig. 5(c)]. A level-set method was applied to obtain a smooth binary retinal mask.⁴⁵ If there was an extra tissue above the retina, the method sometimes created multiple regions that were not connected. Therefore, we classified the region with maximum volume as the retina. Next, vessels were detected as voxels with large Doppler phase shifts inside the retina mask. The threshold was defined as 2.33 times the standard deviation of Doppler phase shift in the static retinal tissue. A morphological operation was then used to exclude small regions and connect nearby regions in the vessel mask. Because spatial interval between two adjacent B-scans was $25\mu \mathrm{m}$, it was difficult to decide if a region with one to two pixels was noise or signal. Therefore, we removed any region with diameter $\le 50\mu \mathrm{m}$ [Fig. 5(e)].

Fig. 5

Detection of veins in Doppler OCT image. (a) Flow chart for vessel detection; (b) Doppler phase shift and reflectance Image on one en face plane. The arrow is the direction of B-scans and A-scans. (c) Binary mask include retina; (d) masked Doppler OCT image using retinal mask; (e) vessel mask in one en face plane after vessel detection; (f) Doppler OCT image in Fig. 5(e) after phase unwrapping; (g) illustration of two types of vein classification rules: vein (class 1) is the vein inside the disc, vein (class 2) is the vein outside the disc. At depth showed in this figure, the en face plane cross with the vein twice; (h) two type of veins in same en face plane; and (i) veins (class 1) detected in one en face plane (blue circle).

Phase wrapping may change the sign of Doppler phase shift in the center of the vessel. Thus, an unsupervised phase unwrapping algorithm^46,47 was applied before the vein and artery classification. The phase unwrapping algorithm was applied to the two-dimensional image in each en face plane separately. The phase unwrapping was only applied to pixels inside the retina mask [Figs. 5(d) and 5(f)].

After the vessels were detected and phase unwrapping algorithm was applied, veins and artery were classified. From the vessel mask, we excluded arteries and veins which were double counted in the same en face plane. Veins were identified as vessels with flow directing to the optic disc in the en face view. We illustrated the vessel type decision in Fig. 5(g). A typical vein had three segments D1 to D3 and the corresponding Doppler phase shifts were negative, close to 0, and positive. The algorithm in the previous paragraph detected two vessel segments: the first one corresponds to D1, or class 1, where the Doppler phase shift slope of vessel segment is negative; and the other one corresponds to D3, or class 2, where Doppler phase shift and slope of vessel segment is positive [Fig. 5(h)]. Vessel segment D2 cannot be detected because the Doppler phase shift is too weak. In arteries, vessel segment D1 had positive Doppler phase shift and negative slope of vessel segment, whereas vessel segment D3 has negative Doppler phase shift and positive slope of vessel segment. We only considered vessel segments with both negative Doppler phase shift and negative slope of vessel segment as valid venous masks for TRBF measurement [Fig. 5(i)].

2.7.

Multiple Planes En Face Doppler Method

Due to artifacts in spectral OCT (such as fringe washout, phase wrapping, or weaker reflectance in deeper tissue), it was difficult to find one artifact-free plane that includes all veins. To solve this problem, we determined an optimized en face plane for each vein branch and then summated the flow of all vein branches (Fig. 6).

Fig. 6

Illustration of multiple plane en face Doppler method. (a)—(c) Three en face planes were selected for comparison. Each vein had an optimized en face plane in which the Doppler phase shift was strong enough but without any artifact. Blue circle is the valid vein mask. Red circle means the particular vessel obtained higher flow measurement than other planes. (d) When two branches merged to one vessel in deeper en face plane, the summation of branches [indicated by the upper two circles in Fig. 6(a)] will be compared to the flow in the merged vessel [indicated by the blue circle in Fig. 6(d)].

We obtained the maximum blood flow on the optimized en face plane for a vein. On the other en face plane, artifacts will reduce the integration of Doppler phase shift in the vessel area and then create lower measurement than the optimized en face plane [Figs. 6(a)—6(c)]. TRBF was summated from the optimized value of each vein. If two vein branches merged to a root vein in deeper level, the summation of branches was compared to the flow in the root vein. Either root or branch veins are selected in the later depending on which one is larger [Fig. 6(d)]. To reduce the effect of cardinal pulsation, the TRBF is also averaged from all volumes in one scan.

The pixel size of the en face plane is proportional to the axial length of the eye. Since axial length is variable, a compensation factor is necessary to correct TRBF using a default pixel size. The compensation factor $c={(L/L_d)}^2$, where $L$ is the axial length of the test eye and $L_d$ is the default axial length. We used $L_d=24\mathrm{mm}$ in this study. This compensation method for vessel area was consistent to methods used in other imaging modalities.⁴⁸

2.8.

Participants

To test the proposed TRBF measuring algorithm, participants were selected from functional and structural OCT for glaucoma study in Casey Eye Institute. The research protocol was approved by the institutional review boards at Oregon Health and Science University and carried out in accordance with the tenets of the Declaration of Helsinki. Written informed consent was obtained from each subject.

In this study, normal participants had intraocular pressures of $<21\mathrm{mm}$ Hg for both eyes, a normal Humphrey Swedish Interactive Threshold Algorithm (SITA) 24-2 standard visual field with mean deviation (MD) and pattern standard deviation within 95% limits of the normal reference. They also had a glaucoma hemifield test within 97% limits, a central corneal thickness $\ge 500\mu \mathrm{m}$, a normal appearing ONH, a normal nerve fiber layer, an open anterior chamber angle as observed by gonioscopy, and no history of chronic ocular or systemic corticosteroid use. The glaucoma participants had at least one eye which had ONH or nerve fiber layer defects visible on slit-lamp biomicroscopy defined as one of following: diffused or localized rim thinning, disc (splinter) hemorrhage, vertical cup/disc ratio greater than the fellow eye by $>0.2$, or a notch in the rim detected on baseline dilated fundus examination and confirmed by masked reading of stereo disc photographs. People were not enrolled if they had other eye diseases, previous intraocular surgery, age $<40$ or $>79$ years, or best corrected visual acuity worse than 20/40.