Doppler optical coherence tomography (OCT) can noninvasively image the flow velocity component along the light beam direction and reveal cross-sectional flow distribution in individual vessels.¹ It enables the calculation of absolute retinal blood flow by either dual-circular² or en face Doppler scanning protocols.³ The dual-circular scan samples two concentric circles near the optic disc and transects all major vessels simultaneously to obtain blood flow.² The en face Doppler scan³ consists of a volumetric scan over the optic disc region and integrates blood flow velocity across the en face plane. This avoids the necessity of vessel Doppler angle determination but requires many more sampling points (A-lines).

In Doppler OCT, blood flow velocity can be calculated from the phase shift between adjacent A-lines or B-scans.¹ To image the major retinal arteries or veins around the optic disc in rodents or humans, the Doppler signal is usually generated by adjacent A-lines, according to the currently available system speeds. Since the phase shift is mathematically restricted to $[-\pi ,\pi ]$,⁴ phase wrapping frequently occurs and confounds the determination of blood flow velocity. Although reducing the time interval of adjacent A-lines or adjusting the Doppler angle closer to 90 deg could enlarge the maximum detectable flow speed ${\nu }_m$ to avoid phase wrapping, these operations will bring additional problems, such as insufficient reflectance and difficult alignment.²

Phase unwrapping is a technique that intentionally expands the phase range to beyond $[-\pi ,\pi ]$. Due to the severe noise level in Doppler OCT, this procedure was initially performed by manually inspecting each frame and isolating the region affected by phase wrapping.⁵ Later, Wang et al.⁶ proposed an automated method that takes advantage of the continuity of blood flow velocity along the A-line profile to correct the phase wrapping. However, correction failure on vessels and miscorrection on tissue still occasionally occur.

In this letter, we propose an automated technique based on the magnitude of the phase-shift gradient to correct phase wrapping. The algorithm works under the presumption that the blood flow velocity within a blood vessel has a parabolic profile, with the highest velocity at the center and lowest velocity near the walls.^6,7 Thus, the phase-wrapping region in cross-sectional Doppler OCT image tends to be limited to a closed area within the vessel. The magnitudes of phase-shift gradients at the edge of this region are extraordinarily large due to the reverse of phase direction, providing an indicator of the presence of phase wrapping, which may be useful for correcting this problem. The effectiveness of this algorithm is demonstrated here by dual-circular Doppler OCT in the rodent retina and en face Doppler OCT in human retina.

Rodent retinal images were acquired by a custom-built visible-light OCT system with a 50-kHz sampling rate, center wavelength at 565 nm, and axial resolution of $1.2\mu \mathrm{m}$ in tissue.⁸ Brown Norway rats were anesthetized with 2.5% isoflurane and the pupil dilated with 0.5% tropicamide ophthalmic solution. All experimental procedures were approved by the Institutional Review Board/Ethics Committee as well as the Institutional Animal Care and Use Committee (IACUC) of the Oregon Health and Science University.

After locating the optic disc [Fig. 1(a)], an ultradense dual-circle (inner radius $r_1=0.6\mathrm{mm}$ and outer radius $r_2=0.8\mathrm{mm}$) scanning protocol was performed, consisting of 4096 A-lines [Fig. 1(b)] around the optic disc for each circle with 30 repeated frames at each location. The scanning densities were $0.9\mu \mathrm{m}/\text{line}$ for the inner circle and $1.2\mu \mathrm{m}/\text{line}$ for the outer circle.

Fig. 1

(a) Mean projection en face structural image of a rat retina, acquired with the prototype vis-OCT system. Dual-circle scanning positions were marked with white dashed circles, with arrows indicating the starting points of scanning ($\theta =0$). Vessels were numbered starting from $\theta =0$ along the counterclockwise scanning direction. The field of view is $2.2\times 2.2{\mathrm{mm}}^2$. $\text{Scale bar}=200\mu \mathrm{m}$. (b) Averaged cross-sectional structural images of the inner circular scan. The image covers $500\mu \mathrm{m}$ in the depth direction and 360 deg in the $\theta$ direction, covering all major retinal vessels. The shadow artifacts underneath superficial vessels can be used to identify the approximate vessel cross section (yellow dash box) for phase unwrapping. (c) A single-frame Doppler phase-shift image of the inner circular scan. Note that phase wrappings exist in some vessels $({\mathrm{V}}_1,{\mathrm{V}}_2)$, as evidenced by the phase jumping ($\text{red}\to \text{blue}$ or $\text{blue}\to \text{red}$) inside the vessel cross-sectional regions.

To reveal the structural images [Fig. 1(b)], interferograms were successively processed by $\lambda -k$ interpolation, DC removal, Hilbert transform, dispersion compensation, and fast Fourier transform. The initial Doppler phase shift was obtained by calculating the angle between adjacent A-lines (complex signal). After that, bulk phase motion between A-lines was removed by standard deviation methods⁹ to obtain the phase shift relevant to in situ blood flow only [Fig. 1(c)].

The horizontal position of each major vessel ($\theta$ direction) was identified by the absorption shadow artifact underneath [Fig. 1(b)] and their depth positions ($z$ direction) were determined by the center of mass of the averaged central vessel A-lines in each B-scan. A $150\text{-}\mu \mathrm{m}\times 13$-deg rectangular area centered at determined coordinates was then used to cover the entire cross section of each vessel [yellow boxes in Fig. 1(b)].

A $5\times 5$ median filter was applied to the regions of interest to remove noise due to abrupt changes in values (${\phi }_0$ in Fig. 2). After further smoothing by a two-dimensional (2-D) Gaussian filter ($\sigma =3$), the magnitude of the phase-shift gradient $G_{\text{mag}}$ [Eq. (1), Fig. 2] by combining the gradients (“Sobel” method) in the $\theta$ direction $G_{\theta }$ and $z$ direction $G_z$ was calculated as

Fig. 2

Phase unwrapping of a vein (${\mathrm{V}}_1$) and an artery (${\mathrm{V}}_2$) affected by phase wrapping as well as a vein (${\mathrm{V}}_9$) unaffected by phase wrapping. First column: phase shift ${\phi }_0$ after median filtering shows that, in the presence of phase wrapping, the phase jumps from negative to positive (blue to red) inside the vein (${\mathrm{V}}_1$) and from positive to negative (red to blue) inside the artery (${\mathrm{V}}_2$). Second column: phase gradient magnitude $G_{\text{mag}}$ was extraordinarily large at the edge of phase wrapping, compared to the vessel periphery and the vessel (${\mathrm{V}}_9$) that did not have phase wrapping. Third column: binary image BW of the phase-wrapping regions. Note their presence in the vessels with phase wrapping and their absence in the vessel without phase wrapping. Fourth column: phase-unwrapped Doppler signal ${\phi }_{puw}$ as determined by the proposed algorithm. The phase shift was corrected and shown within the range between $-2\pi$ and $2\pi$. Fifth column: parabolic flow velocity profile retrieved from unwrapped phase shift after registering and averaging the results from all repeated frames. Image size: $150\mu \mathrm{m}(\mathrm{z})\times 13\mathrm{deg}(\theta )$.

At the edge of the phase-wrapping region, the phase shift reverses from $-\pi$ to $\pi$ or $\pi$ to $-\pi$, causing a gradient magnitude ($\sim 2\pi$) significantly larger than that suggested by the regular velocity profile. However, the median and Gaussian filters degrade the gradient magnitude. Therefore, the superthreshold was set to $\pi /2$ to obtain the edge of the phase-wrapping region (bright contour in $G_{\text{mag}}$ in Fig. 2). Ideally, the edge contour should be a closed curve due to the parabolic flow velocity in vessels. However, breaks in the contour occur due to phase discontinuity. Morphological closing (an operation that will remove isolated below-threshold pixels within a region) was applied to acquire a closed phase-wrapping region (BW in Fig. 2).

The phase correction factor ($2\pi$ or $-2\pi$) was determined by the phase shift outside the phase-wrapping region. This ring region was easily generated by dilating the phase-wrapping region BW and then subtracting the original BW from the dilated region. An averaged phase ${\phi }_p$ was calculated over the annular region along the vessel wall as an indicator of flow direction. The unwrapped phase signal ${\phi }_{puw}$ was then obtained by adding $2\pi$ if ${\phi }_p$ was positive (${\mathrm{V}}_2$ in Fig. 2) or subtracting $2\pi$ if ${\phi }_p$ was negative (${\mathrm{V}}_1$ in Fig. 2). The full phase unwrapping for each vessel at each frame takes $<0.1\mathrm{s}$ to accomplish in the platform of MATLAB R2018b with Intel i7-7700 CPU processor

After that, the Doppler B-scan frames of each vessel were rigidly registered using the reflectance images by a 2-D cross-correlation method,¹⁰ and the phase shifts were averaged over the registered B-scans to account for the cardiac cycle. The Doppler angle $\alpha$ of each vessel was calculated from the difference of the positions of vessel centers (determined by the center of mass in ${\phi }_{puw}$ images) between the inner and outer circular scans. The parabolic profiles of the blood flow velocities (Fig. 2) were calculated according to $\nu =f_s{\lambda }_0\phi /[4\pi n\cos (\alpha )]$, where the sampling rate $f_s$ is 50 kHz, the center wavelength ${\lambda }_0$ is 565 nm, and refractive index $n$ is 1.4. Mean blood flow velocity $\bar{\nu }$ was calculated by averaging the velocity distribution profiles on the vessel region (acquired by setting a $1\text{-}\mathrm{mm}/\mathrm{s}$ threshold to the velocity distribution profiles), with vessel diameter $H$ in B-scans obtained by the maximum axial span of the vessel area. After that, vessel cross-sectional area $S$ was calculated from the vessel diameter $H$ and the Doppler angle $\alpha$ as $S=\pi {[H\sin (\alpha )/2]}^2$ under the presumption of a cylindrical shape for a vessel. The blood flow $F$ was calculated as the product of mean flow velocities $\overline{\nu }$ over the cross-sectional area $S$ of each vessel as follows:

The measurements for retinal vessels shown earlier are listed in Table 1. For the vessels affected by phase wrapping, the mean phase shifts were higher than that with smaller Doppler angles. The measured total retinal blood flow calculated from venous flow was $8.53\pm 0.87\mu \mathrm{L}/\min$ ($\text{rats}=6$) with mean flow velocities ranging from 7 to $12\mathrm{mm}/\mathrm{s}$. These results are consistent with previously reported values.¹¹

Table 1

Mean phase shift φ, Doppler angle α, mean velocity ν, and blood flow f for the major vessels in the rat retina shown in Fig. 1.

Human retinal images were acquired using a commercial spectral-domain OCT system (RTVue-XR Avanti; Optovue, Inc., Fremont, California) with a 70-kHz scanning rate, center wavelength at 840 nm, and an axial resolution of $5\mu \mathrm{m}$. To achieve en face Doppler measurements, a dense volumetric scan consisting of 500 A-lines/B-scans and 78 B-scans was performed on a $2\times 2{\mathrm{mm}}^2$ region around the optic disc [Fig. 3(a)], followed by 5 repeated volume determinations to compensate for the cardiac cycle. The whole acquisition time is about 3 s.

Fig. 3

Doppler phase unwrapping in a human retina. (a) Human retina centered at optic disc region ($2\times 2{\mathrm{mm}}^2$) by commercial infrared OCT system. (b) Reflectance B-scan image in the position marked by white dash line in (a). (c) Doppler signal ${\phi }_0$, with three vessels marked by red arrows. (d) Bright contours (red arrows) in phase gradient magnitude map $G_{\text{mag}}$ indicating two vessels affected by phase wrapping. (e) Corresponding binary image BW for the phase-wrapping regions. (f) Corresponding phase-unwrapped Doppler signal ${\phi }_{puw}$. (g) Volumetric Doppler phase-shift signal with veins is identified by flow direction.

The blood absorption attenuation in infrared is weaker than that in visible light; however, the much larger diameter in depth in human retinal vessels causes a dark shadow similar to that in vis-OCT [Fig. 3(b)]. Since vessels come almost vertically from the optic disc, the Doppler angles are far away from 90 deg for arteries and veins thus phase wrapping is easily seen in phase-shift images [Fig. 3(c)].

For the en face Doppler OCT, the phase gradient magnitude map $G_{\text{mag}}$ can be calculated from the whole B-scan. Two vessels appear bright in Fig. 3(d), indicating that they are affected by phase wrapping within the regions delimited by BW [Fig. 3(e)]. The phase was unwrapped successively inside the phase-wrapping regions using the algorithm described earlier [Fig. 3(f)].

The differentiation between arteries and veins was assisted by three-dimensional visualization in this scanning protocol [Fig. 3(g)]. Since the Doppler angle of vessels varies dramatically within the optic disc region, the phase shift for one vessel can be changed between negative (blue) and positive (red), with intermediate breaks where phase shift is near zero at a Doppler angle close to 90 deg.

Usually, the total retina blood flow in humans is calculated in veins. This is because, as opposed to arteries, veins (1) have steady flow during systole and diastole and (2) are straight and easy to find a satisfied depth plane. As shown in Fig. 3(g), the superior ophthalmic vein trunk (vein 1), nasal- (vein 2), and temporal- (vein 3) inferior ophthalmic vein branches, as well as two ciliary veins (veins 4 and 5) were identified by their morphology and flow directions. By contrast, the arteries broke into multiple slices and were hard to delineate for flow measurement.

When the velocity is too high, multiple-fold phase wrappings may occur, which can produce multiple nested bright contours in the phase gradient magnitude map. The phase then should be unwrapped by correcting the phase in the outer contour first and then the phase in the inner contour. However, in these cases, fringe washout can be severe and make it difficult to obtain the phase shift from adjacent A-lines. So in practice, multiple-fold phase wrapping of multiple vessels should still be avoided.

Limitations still exist in this algorithm. First, the proposed phase-unwrapping algorithm is strictly specialized for determining blood flow, and thus it is not suitable as a generic algorithm that can be easily adapted to other applications. However, the idea of introducing image processing techniques to Doppler phase unwrapping may be generally beneficial. Second, if the flow profile is nonparabolic (for instance because of a sharp bend in the vessel), errors in the detection of the wrapping boundary might occur. For example, the velocity profile in the middle vessel in Fig. 3(c) is irregular, causing an imperfectly closed phase-wrapping region. If the irregularity were more severe, the phase-unwrapping process would fail.

In summary, under the assumption that phase wrapping predominantly occurs in the vessel center, we used the phase gradient magnitude to indicate the presence of and correct for phase wrapping. The algorithm was verified using a dual-circle protocol with a visible-light prototype in the rat retina and an en face Doppler scanning protocol with an infrared commercial OCT system in the human retina.