## Abstract

Retinal blood flow quantification by retinal vessel segmentation with Doppler optical coherence angiography is presented. Vessel diameter, orientation, and position are determined in an *en face* vessel image and two representative cross-sectional flow images of the vessel. Absolute blood flow velocity is calculated with the help of the measured Doppler frequency shift and determined vessel angle. The volumetric flow rate is obtained with the position and the region of the vessel lumen. The volumetric blood flow rate of retinal arteries before and after a bifurcation is verified in a healthy human eye.

© 2008 Optical Society of America

Retinal blood flow measurement is required to determine the effect of hemodynamic abnormalities in retinal pathologies. There exist several blood flow measurement techniques that are applied to retinal vessels. Dye dilution is a technique that involves tracing fluorescent dyes flowing in the retinal vessel by monitoring the fluorescence intensity. However, the dye injection is not comfortable for patients. Although laser Doppler velocimetry (LDV) can measure the volumetric flow rate of the retinal vessel [1], it requires dilation to collect two scattered beams with different scattering angles.

Doppler optical coherence tomography (OCT) is a technique used to observe cross-sectional flow images, and retinal blood flow distribution has been imaged. By using a spectral-domain detection technique, high-speed cross-sectional retinal blood flow imaging is achieved [2, 3]. Furthermore, the three-dimensional (3D) retinal vasculature can be observed [4, 5].

To obtain the flow quantification with the Doppler OCT, several techniques have been developed, e.g., axial and transversal flow speeds measurement [6] and multiangle detection techniques [7, 8]. However, these techniques have several shortcomings for *in vivo* retinal blood flow measurement. For example, to measure the transversal velocity, a large number of samples is required. The multiangle detection technique may require dilation, as LDV does. In addition, it results in the limited axial range in spectral-domain OCT (SD-OCT). Recently, vessel segmentation techniques that calculate the vessel orientation were presented. Wang *et al.* [9] presented a total blood flow volume measurement with sampling all retinal arteries or veins around the optic nerve head. The vessel orientations are obtained from vectors connecting vessel centers between two successive cross sections. In this technique, a straight portion of the vessel is aimed, and these two cross sections should not be nearly parallel to the vessel. Full 3D processing was performed by Michaely *et al.* [10]. Although this method is comprehensive, it requires a long computational time, since an OCT volume set has a large number of voxels. In this Letter, we demonstrate a retinal blood flow quantification with a two-step computationally inexpensive retinal vessel segmentation with a two-dimensional (2D) *en face* retinal vessel image and 2D cross-sectional flow image.

For the absolute blood flow speed measurement with the Doppler OCT, the Doppler angle is required, and for the volumetric flow measurement, vessel lumen determination is needed. For that purpose, vessel segmentation is applied to determine the vessel parameters, i.e., the vessel orientation and diameter.

Since the retinal vasculature is distributed nearly parallel to the retinal surface, *en face* retinal vessel images illustrate the representative retinal vasculature. Thus, the vessel segmentation was performed in two steps. The 3D geometric parameters of the vessel are indicated in Fig. 1 . The first step is the vessel azimuth *θ*, diameter *D*, and lateral position detection in the *en face* projection image *E*. The second step is the vessel zenith *ϕ* and axial position detection in the cross-sectional vessel images (planes A and P in Fig. 1). These vessel parameters are then used to determine the region of the vessel lumen *S*. This algorithm does not require processing for the entire 3D volume set; thus, the computational time is shortened.

The second derivatives of the *en face* vessel image are used to determine the parameters of the curvilinear structure, namely, retinal vessels. For the segmentation of several sizes of vessels, a multiscale technique [11] is used. The filtering operator is defined with second derivatives of a Gaussian function with a standard deviation *σ*, as $L={[\partial \u2215\partial x,\partial \u2215\partial y]}^{T}[\partial \u2215\partial x,\partial \u2215\partial y]G(\mathbf{r},\sigma )$. Applying this line filter to the vessel image *I*, we have a $2\times 2$ matrix of the second derivatives of the image $\mathbf{H}(x,y)=L\otimes I(x,y)$. Vesselness *F* is calculated from the eigenvalues of the matrix **H**, ${\Lambda}_{1}$ and ${\Lambda}_{2}$
$({\Lambda}_{1}>{\Lambda}_{2})$, as follows: $F=\sqrt{\mid {\Lambda}_{2}\mid (\mid {\Lambda}_{2}\mid +{\Lambda}_{1})}\phantom{\rule{0.2em}{0ex}}({\Lambda}_{2}<{\Lambda}_{1}\u2a7d0)$, $\sqrt{\mid {\Lambda}_{2}\mid (\mid {\Lambda}_{2}\mid -\alpha {\Lambda}_{1})}\phantom{\rule{0.2em}{0ex}}({\Lambda}_{2}<0<{\Lambda}_{1}<\mid {\Lambda}_{2}\mid \u2215\alpha )$, or 0 (otherwise). These conditions are determined to discriminate a bloblike shape $(\mid {\Lambda}_{1}\mid \sim \mid {\Lambda}_{2}\mid \u2aa20)$ from a curvilinear structure and make it sensitive to the discontinuity of the line structure $({\Lambda}_{2}<0<{\Lambda}_{1})$. The value of *α* that determines the sensitivity to discontinuity is set to 0.25. Here, the vessels in the image are considered to possess a brighter intensity than that of the surroundings. Usually, retinal vessels cast shadows in an *en face* projection image of an OCT volume set, so that an inverted *en face* image is used. The eigenvectors indicate the orientation of the line structure, namely, azimuth angle of the vessel *θ*. For a multiscale detection, this line filter is applied with several different deviations of the Gaussian function. Vesselness is then normalized with the deviation of the Gaussian function *σ* as $\sigma F\left(\sigma \right)$. The vesselness at the applied point is defined as the maximum normalized vesselness ${F}^{\prime}=\mathrm{max}\left[\sigma F\left(\sigma \right)\right]$, and the vessel diameter is defined as $D=4{\sigma}_{\mathrm{max}}$, where ${F}^{\prime}={\sigma}_{\mathrm{max}}F\left({\sigma}_{\mathrm{max}}\right)$. The measured vessels are selected semiautomatically. The retinal vessels in the *en face* image are determined with a top-hat filter, followed by thresholding. The operator indicates a line cross over a retinal vessel. Then, the point along that line having the maximum vesselness is defined as the center of the vessel.

According to the azimuth angle *θ* and the central position of the vessel, a cross-sectional flow image (plane A in Fig. 1) line intersecting the vessel centerline is created from the 3D flow volume set. Line filtering with the standard deviation of the Gaussian function ${\sigma}_{\mathrm{max}}$ is applied to the unsigned flow image. The zenith angle $\left(\varphi \right)$ of the vessel corresponds to the line filter output at the center of the vessel. To determine the vessel position and the region of the vessel lumen, the depth of the vessel should be known. A cross-sectional bidirectional flow image intersecting the vessel (plane P in Fig. 1) is created from the 3D flow volume set. The correlation between this cross section and a circle with a diameter *D* is applied. The position of the maximum correlation point indicates the depth of the vessel. The zenith angle of the vessel can be considered as the Doppler angle when all the blood flow vectors are parallel.

With the help of the Doppler angle *ϕ* and Doppler frequency shift $\Delta f$, the blood flow velocity can be corrected, except in the case of $\varphi \approx 90\xb0$. Although this exceptional Doppler angle must be avoided in practice, it does not frequently occur. The absolute blood flow velocity *V* is obtained as $V={\lambda}_{0}\Delta f\u2215\left(2n\phantom{\rule{0.2em}{0ex}}\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\varphi \right)$, where ${\lambda}_{0}$ is the central wavelength of the light source and *n* indicates the refractive index of the sample. According to the obtained vessel parameters, an elliptical vessel lumen *S* in the cross-sectional flow image (plane P in Fig. 1) is estimated. The blood flow rate is obtained by the integration of the lateral blood flow velocity component, which is calculated from the Doppler frequency shift and the Doppler angle within the vessel lumen *S* as $J={\sum}_{S}{\lambda}_{0}\Delta f\u22152n\phantom{\rule{0.2em}{0ex}}\mathrm{tan}\phantom{\rule{0.2em}{0ex}}\varphi \Delta s$, where $\Delta s$ is the area of one pixel in the cross-sectional flow image.

For the verification of this method, blood mass conservation is confirmed at a bifurcation of a retinal artery. The right eye of a healthy volunteer was scanned by an ultrahigh-resolution SD-OCT system [12] with $27.7\text{\hspace{0.17em}}\mathrm{kHz}$ of axial scan rate and $\sim 3\text{\hspace{0.17em}}\mu \mathrm{m}$ of axial resolution in the tissue. A 3D volume set with $196\times 28$ axial scans (corresponding to $1\text{\hspace{0.17em}}\mathrm{mm}\times 1\text{\hspace{0.17em}}\mathrm{mm}$) was obtained within $225\text{\hspace{0.17em}}\mathrm{ms}$. Since the voxel size is anisotropic, the segmentation accuracy and resulting measurement of blood flow will depend on the azimuth angle. Twenty volumes are sequentially scanned so that the total measurement time is $4.5\text{\hspace{0.17em}}\mathrm{s}$. Since almost all of the retinal vessels are nearly perpendicular to the incident beam, the Doppler frequency shifts can be relatively small. This results in a significant error in the vessel segmentation algorithm. To avoid such a situation, an artery, rather than a vein, is used since the Doppler frequency shift of the arteries is larger than that of the veins. For Doppler shift flow imaging, the Kasai correlator is applied to the complex OCT image. To avoid a $2\pi $-phase ambiguity in the flow image, which causes the limited Doppler frequency range, a 2D-phase unwrapping algorithm [13] is applied to a portion of the flow image wherein the phase wrapping has occurred. Motion artifacts of the flow images are removed by a histogram-based method, and image distortion due to the bulk motion is numerically corrected [4]. To obtain a high-contrast *en face* retinal vessel image, the maximum log-intensity projection of the 3D retinal OCT volume set just around the retinal pigment epithelium (RPE) is calculated (Fig. 2A ).

The multiscale line filtering results of one measurement are shown in Figs. 2B–2D. Accordingly, the vessel segmentation is performed and results are shown in Fig. 3 . The computational time of the vessel segmentation and flow calculation for each volume (three vessels) are $\sim 3\text{\hspace{0.17em}}\mathrm{min}$ with a personal computer (AMD Athlon 64 $3500+$, $2\text{\hspace{0.17em}}\mathrm{GB}$ RAM). A semiautomatic registration of the *en face* image and a fundus photograph (Fig. 3A) is performed, and the area corresponding to the scanning region is indicated as the white box. The red lines in Fig. 3B indicate the cross-sectional planes parallel to the retinal vessels, and the blue lines are the perpendicular cross sections. At the vessel 2, the azimuth angle has the small tilt, although the vessel seems to be nearly horizontal. However, the corresponding fundus photography (Fig. 3A) shows the slight azimuth angle of vessel 2. It can be considered that the multiscale line filtering is sensitive to the subpixel ordered structure. The cross-sectional flow images parallel and perpendicular to the retinal vessel are shown in Figs. 3C–3H. These images consist of the nearest voxels from the 3D flow volume set according to *en face* positions. The determined vessel boundaries are indicated with ellipses in Figs. 3C, 3E, and 3G. The red lines in Figs. 3D, 3F, and 3H are the tangential lines to the centerline of the vessels at the center of the images.

The blood flow parameters are obtained for all of the 20 volume sets and the mean of the peak velocity ${V}_{\text{peak}}$, diameter *D*, and volumetric flow rate *J* are shown in Table 1 . The signs of the blood flow velocity and rate indicate the direction of flow, incoming to (+) or outgoing from (−) the bifurcation. The plot of the volumetric flow rate of vessel 1 versus the total flow rate of vessels 2 and 3 is shown in Fig. 4 . The measured incoming and outgoing blood volumes are in direct proportion. The intraclass correlation coefficient is 0.74. The mean of the total outgoing blood flow rate, $3.41\text{\hspace{0.17em}}\mu \mathrm{l}\u2215\mathrm{min}$, is comparable to the incoming blood flow rate, $3.23\text{\hspace{0.17em}}\mu \mathrm{l}\u2215\mathrm{min}$. Since the blood flow in the arteries are measured, the standard deviations of the blood flow rates are relatively large as compared to the mean flow rates due to pulsation.

In conclusion, the vessel parameters—orientation and diameter—are determined with a two-step 2D image-based processing method. The retinal blood flow velocity and volumetric flow rate are quantitatively obtained. To the best of our knowledge, this is the first study of conservation of blood flow at a bifurcation of the retinal artery with a cross-sectional flow imaging technique.

The authors acknowledge the helpful discussion with Masahiro Miura, Department of Ophthalmology, Tokyo Medical University. This research is supported by a grant-in-aid for scientific research 18360029 and the Japan Society for the Promotion of Science (JSPS) fellows 18⋅3827 from JSPS, Japan Science and Technology Agency, and the Special Research Project of Nanoscience at the University of Tsukuba. Shuichi Makita is a JSPS research fellow.

**1. **C. E. Riva, J. E. Grunwald, S. H. Sinclair, and B. L. Petrig, Invest. Ophthalmol. Visual Sci. **26**, 1124 (1985).

**2. **R. Leitgeb, L. Schmetterer, W. Drexler, A. Fercher, R. Zawadzki, and T. Bajraszewski, Opt. Express **11**, 3116 (2003). [CrossRef] [PubMed]

**3. **B. R. White, M. C. Pierce, N. Nassif, B. Cense, B. H. Park, G. J. Tearney, B. E. Bouma, T. C. Chen, and J. F. de Boer, Opt. Express **11**, 3490 (2003). [CrossRef] [PubMed]

**4. **S. Makita, Y. Hong, M. Yamanari, T. Yatagai, and Y. Yasuno, Opt. Express **14**, 7821 (2006). [CrossRef] [PubMed]

**5. **A. H. Bachmann, M. L. Villiger, C. Blatter, T. Lasser, and R. A. Leitgeb, Opt. Express **15**, 408 (2007). [CrossRef] [PubMed]

**6. **D. Piao and Q. Zhu, Appl. Opt. **42**, 5158 (2003). [CrossRef] [PubMed]

**7. **C. J. Pedersen, D. Huang, M. A. Shure, and A. M. Rollins, Opt. Lett. **32**, 506 (2007). [CrossRef] [PubMed]

**8. **Y. C. Ahn, W. Jung, and Z. Chen, Opt. Lett. **32**, 1587 (2007). [CrossRef] [PubMed]

**9. **Y. Wang, B. A. Bower, J. A. Izatt, O. Tan, and D. Huang, J. Biomed. Opt. **12**, 041215 (2007). [CrossRef] [PubMed]

**10. **R. Michaely, A. H. Bachmann, M. L. Villiger, C. Blatter, T. Lasser, and R. A. Leitgeb, J. Biomed. Opt. **12**, 041213 (2007). [CrossRef] [PubMed]

**11. **T. Koller, G. Gerig, G. Szekely, and D. Dettwiler, in Proceedings of the Fifth International Conference on Computer Vision (IEEE, 1995), p. 864.

**12. **Y. Hong, S. Makita, M. Yamanari, M. Miura, S. Kim, T. Yatagai, and Y. Yasuno, Opt. Express **15**, 7538 (2007). [CrossRef] [PubMed]

**13. **D. Ghiglia and L. Romero, J. Opt. Soc. Am. A **11**, 107 (1994). [CrossRef]