|
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.1 It enables the calculation of absolute retinal blood flow by either dual-circular2 or en face Doppler scanning protocols.3 The dual-circular scan samples two concentric circles near the optic disc and transects all major vessels simultaneously to obtain blood flow.2 The en face Doppler scan3 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.1 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 ,4 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 to avoid phase wrapping, these operations will bring additional problems, such as insufficient reflectance and difficult alignment.2 Phase unwrapping is a technique that intentionally expands the phase range to beyond . 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.5 Later, Wang et al.6 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 in tissue.8 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 and outer radius ) 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 for the inner circle and for the outer circle. To reveal the structural images [Fig. 1(b)], interferograms were successively processed by 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 methods9 to obtain the phase shift relevant to in situ blood flow only [Fig. 1(c)]. The horizontal position of each major vessel ( direction) was identified by the absorption shadow artifact underneath [Fig. 1(b)] and their depth positions ( direction) were determined by the center of mass of the averaged central vessel A-lines in each B-scan. A -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 median filter was applied to the regions of interest to remove noise due to abrupt changes in values ( in Fig. 2). After further smoothing by a two-dimensional (2-D) Gaussian filter (), the magnitude of the phase-shift gradient [Eq. (1), Fig. 2] by combining the gradients (“Sobel” method) in the direction and direction was calculated as At the edge of the phase-wrapping region, the phase shift reverses from to or to , causing a gradient magnitude () 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 to obtain the edge of the phase-wrapping region (bright contour in 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 ( or ) 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 was calculated over the annular region along the vessel wall as an indicator of flow direction. The unwrapped phase signal was then obtained by adding if was positive ( in Fig. 2) or subtracting if was negative ( in Fig. 2). The full phase unwrapping for each vessel at each frame takes 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,10 and the phase shifts were averaged over the registered B-scans to account for the cardiac cycle. The Doppler angle of each vessel was calculated from the difference of the positions of vessel centers (determined by the center of mass in images) between the inner and outer circular scans. The parabolic profiles of the blood flow velocities (Fig. 2) were calculated according to , where the sampling rate is 50 kHz, the center wavelength is 565 nm, and refractive index is 1.4. Mean blood flow velocity was calculated by averaging the velocity distribution profiles on the vessel region (acquired by setting a threshold to the velocity distribution profiles), with vessel diameter in B-scans obtained by the maximum axial span of the vessel area. After that, vessel cross-sectional area was calculated from the vessel diameter and the Doppler angle as under the presumption of a cylindrical shape for a vessel. The blood flow was calculated as the product of mean flow velocities over the cross-sectional area 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 () with mean flow velocities ranging from 7 to . These results are consistent with previously reported values.11 Table 1Mean 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 . 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 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. 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 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. DisclosuresAcner Camino and Yali Jia have a significant financial interest in Optovue Inc. These potential conflicts of interest have been reviewed and managed by the Oregon Health and Science University. Other authors declare that there are no conflicts of interest related to this letter. AcknowledgmentsThis work was supported by Grant Nos. R01 EY027833, DP3 DK104397, R01 EY024544, R01 EY010145, and P30 EY010572 from the National Institutes of Health, Bethesda, Maryland, and an unrestricted departmental funding grant and William and Mary Greve Special Scholar Award from the Research to Prevent Blindness, New York, New York. ReferencesR. A. Leitgeb et al.,
“Doppler optical coherence tomography,”
Prog. Retinal Eye Res., 41 26
–43
(2014). https://doi.org/10.1016/j.preteyeres.2014.03.004 Google Scholar
Y. Wang et al.,
“In vivo total retinal blood flow measurement by Fourier domain Doppler optical coherence tomography,”
J. Biomed. Opt., 12 041215
(2007). https://doi.org/10.1117/1.2772871 JBOPFO 1083-3668 Google Scholar
O. Tan et al.,
“En face Doppler total retinal blood flow measurement with 70 kHz spectral optical coherence tomography,”
J. Biomed. Opt., 20 066004
(2015). https://doi.org/10.1117/1.JBO.20.6.066004 JBOPFO 1083-3668 Google Scholar
K. Itoh,
“Analysis of the phase unwrapping algorithm,”
Appl. Opt., 21 2470
(1982). https://doi.org/10.1364/AO.21.002470 APOPAI 0003-6935 Google Scholar
B. R. White et al.,
“In vivo dynamic human retinal blood flow imaging using ultra-high-speed spectral domain optical Doppler tomography,”
Opt. Express, 11 3490
–3497
(2003). https://doi.org/10.1364/OE.11.003490 OPEXFF 1094-4087 Google Scholar
Y. Wang et al.,
“Two-dimensional phase unwrapping in Doppler Fourier domain optical coherence tomography,”
Opt. Express, 24 26129
–26145
(2016). https://doi.org/10.1364/OE.24.026129 OPEXFF 1094-4087 Google Scholar
R. A. Leitgeb et al.,
“Real-time assessment of retinal blood flow with ultrafast acquisition by color Doppler Fourier domain optical coherence tomography,”
Opt. Express, 11 3116
–3121
(2003). https://doi.org/10.1364/OE.11.003116 OPEXFF 1094-4087 Google Scholar
S. Pi et al.,
“Angiographic and structural imaging using high axial resolution fiber-based visible-light OCT,”
Biomed. Opt. Express, 8 4595
–4608
(2017). https://doi.org/10.1364/BOE.8.004595 BOEICL 2156-7085 Google Scholar
X. Wei et al.,
“Fast and robust standard-deviation-based method for bulk motion compensation in phase-based functional OCT,”
Opt. Lett., 43 2204
–2207
(2018). https://doi.org/10.1364/OL.43.002204 OPLEDP 0146-9592 Google Scholar
M. Guizar-Sicairos, S. T. Thurman and J. R. Fienup,
“Efficient subpixel image registration algorithms,”
Opt. Lett., 33 156
–158
(2008). https://doi.org/10.1364/OL.33.000156 OPLEDP 0146-9592 Google Scholar
W. Song et al.,
“A combined method to quantify the retinal metabolic rate of oxygen using photoacoustic ophthalmoscopy and optical coherence tomography,”
Sci. Rep., 4 6525
(2014). https://doi.org/10.1038/srep06525 SRCEC3 2045-2322 Google Scholar
|