1.
Introduction
Most retinal diseases are closely related to vasculopathy, such as diabetic retinopathy, retinal vascular occlusive disease, and choroidal neovascularization. Fundus fluorescein angiography (FFA) is widely used in clinical practice and is the gold standard for vascular imaging of the retina and choroid. However, it is an invasive procedure with an associated risk of complications which limits clinical applications.1 In addition, FFA can only provide two-dimensional (2-D) images of the fundus, and the deeper capillary network may not be visualized well by FFA.2
Optical coherence tomography (OCT) is a noninvasive, high-resolution biomedical imaging technology that can provide three-dimensional (3-D) images of the fundus. Doppler OCT (D-OCT) is a functional extension of OCT which can image not only structure but also blood flow.3–5 Compared to the absolute value of flow velocity, the microvascular networks in the retina are also important for diagnosing retinal diseases. A number of methods have been developed to visualize the microvascular networks. A phase-resolved Doppler variance method was first used to map vessels in human skin and brain.6,7 Since then, it has been used in imaging retinal flow.8–10 Yasuno et al.11 and Wang et al.12 used a modified Hilbert transform to achieve high-resolution images of blood flow. Several extensions of D-OCT based on amplitude variance have also demonstrated capabilities of microvascular imaging which are not sensitive to phase noise artifacts, including speckle variance,13,14 correlation mapping,15 split-spectrum amplitude-decorrelation angiography,16 intensity-based Doppler variance (IBDV).9,17,18 Compared with phase-resolved methods, amplitude variance methods do not depend on phase stability. The IBDV method has been demonstrated in a phase instable situation, and the human choroidal blood vessel network was also achieved.9
Compared with FFA which only shows 2-D images of the fundus, OCT can provide the 3-D structure of the fundus. The vascular morphology of intraretinal layers can be obtained by combining OCT angiography and 3-D segmentation of intraretinal layers. It may allow earlier diagnosis and more precise monitoring of some retinal diseases.19,20 There are several approaches that have successfully segmented intraretinal layers in 2-D images, such as the active contour model,21 Dijkstra algorithm,22 and dynamic programming.23 The 3-D surfaces of the retina can be obtained by applying these 2-D based algorithms independently on each image. Compared with 2-D-based methods, 3-D-based segmentation methods can make full use of the information from neighboring B-scans and help improve the accuracy and robustness of the algorithm. A graph cut-based segmentation method was proposed for simultaneous segmentation of multiple 3-D surfaces in the macula24 and optic nerve head.25 However, the minimum cut graph algorithm is computationally expensive which usually requires minutes to complete a 3-D data segmentation. Cheng and Lin26 proposed a fast 3-D dynamic programming expansion method for vessel boundary detection on magnetic resonance imaging sequences. Efficiency was significantly improved with good robustness, and this method can be used in retinal boundary detection.
In this paper, we combine the 3-D dynamic programming method with intensity-based Doppler variance based on swept-source OCT (SS-OCT). The microvascular morphology of intraretinal layers is visualized in vivo in a normal subject by using this method. The microvascular network with the 3-D structure of the retina has potential for earlier diagnosis and precise monitoring in retinal vascular diseases.
2.
Methods
In this study, a 1050-nm SS-OCT system was used to acquire the 3-D data of the fundus. A commercially available swept-source laser (Axsun Technologies Inc., Billerica, Massachusetts) with a center wavelength of 1050 and 100 nm tuning range was used. The configuration of the system has already been described in previous reports.27 The axial resolution was , and the image range was 3.0 mm at the retina (). The interference signal was digitized by a 12 bit 500 MHz data acquisition board, ATS9350 (Alazar Technologies Inc., Pointe-Claire, QC, Canada), and the signal acquisition was triggered by an optical clock generated by the Axsun laser. The signal-to-noise ratio (SNR) of the system was 95 dB with 1.8 mW of sample arm power and a 100-kHz acquisition rate. A 3-dB sensitivity roll-off from 0.1-mm imaging depth to 1.5-mm imaging depth was measured.
Three right eyes (two normal subjects and one high myopia patient) were used to demonstrate the feasibility of the method under the protocols approved by the Institutional Review Board. Each set of 3-D measurements was composed of 200 slices of 400 A-lines over an area of . Each slice was composed of eight sequential B-scans at the same position, which can achieve a high time difference () to improve the sensitivity of the IBDV method. The total acquisition time of a 3-D volumetric data was approximately 6.8 s including the flyback time. The subject was asked to fixate on an internal target during the acquisition.
Custom-built OCT data acquisition and processing software was developed with Microsoft Visual C++ 2012. In order to achieve real-time preview of the acquired data, the graphics processing unit (GPU)-accelerated computing technique was used with NVIDIA’s compute unified device architecture (CUDA, 6.5) toolkit. The SS-OCT preprocessing, including fixed pattern noise removal, spectral shaping, zero-padding, and numerical dispersion compensation, was accelerated by the GPU. Figure 1 shows the flow diagram of the IBDV method combined with 3-D segmentation after obtaining the complex OCT data by preprocessing the interference signal.
Fig. 1
Flow diagram of intensity-based Doppler variance imaging process with three-dimensional (3-D) segmentation.

The intensity-based Doppler variance method9 was used to visualize the retinal vasculature. The time difference () was increased by using the interframe method. Although the motion artifacts were significantly attenuated because of the improvement of acquisition speed, the motion between the B-scans still cannot be ignored, especially the axis motion.28 A subpixel registration algorithm29 was used to alignment sequential B-scans in both the fast transverse direction, , and axial direction, . The variance value () between sequential B-scans in the same position is given by9,18
where is the number of the depth points that are averaged and is the number of the B-scans at the same position. The SNR can be increased by choosing a larger and . In this application, the is the set at 2 and the is the set at 8. denotes the amplitude value which is extracted from complex OCT data. Finally, a threshold based on the histogram analysis of the averaged intensity images is used to remove the noise.In order to obtain the vascular morphology of the intraretinal layers, the boundaries of each layer are required. We used a 3-D expansion of the dynamic programming method26 for automated boundaries detection. The basic idea of the method is to transform the boundary detection problem to an optimization problem that searches for an optimal path. In 2-D images, each pixel is treated as a node in the graph, and the links connecting the nodes are called edges. The values of the links are assigned based on the intensity of the pixels. The start node is a virtual node which connects to the nodes in the first column in the image and the weights are assigned to zero. Correspondingly, the end node is also a virtual node which connects to the nodes in last column. The boundaries can be found by searching the minimum cost of the path from start node to end node. For 3-D expansion, the surfaces of the layers were achieved by sequentially performing the search process in the fast scan plane and slow scan plane. The dynamic programming provides an efficient optimization method to find the minimum cost of the surface.
Consider a volumetric image of size , where the and are the numbers of fast scan and slow scan, and denotes the depth. The iterative cost function in two directions is given by26
for and , where and are the accumulation costs in the plane and the plane, respectively. The initialization of and are assigned by . The parameters and constrain the searching range which can be used to control the smoothness of the surfaces. The parameters and are also used to control the smoothness of the surfaces. For retinal segmentation, the parameters and are set to zero. The volumetric image is the intensity gradient image along the vertical direction which is obtained by the dark-to-light or light-to-dark intensity transition. denotes the weights from node to node and the weights can be calculated as follows22 where is the gradient value of the node which is normalized to values between 0 and 1. The is the minimum weight for system stabilization, which is set at .Figure 1 outlines the core steps of the 3-D segmentation method for the retina. The averaged images were obtained from eight consecutive B-scans after registration. Before segmentation, a 3-D median filter is applied to the volumetric data using a window size of voxels. The initial location of the retina is determined by the inner limiting membrance (ILM) layer and retinal pigment epithelium (RPE) layer. The ILM layer is defined as the first high reflective increase from the inner side of the retina, and the RPE layer can be determined as the highest peak after the ILM layer. Therefore, the initial location of the ILM and RPE layer can be determined from the A-scan profile. After smoothing by using the polynomial method, the search area for the ILM layer can be reduced based on the rudimentary ILM boundary. Finally, the ILM layer is refined by finding the optimal path in a limited region. After obtaining the location of the retina, the limited search region can be determined based on the characteristics of the retinal structure. The OPL/ONL layer and IS/OS layer are the most obvious layers besides the ILM layer, and the intensity transition along the axis direction between two layers are in opposite directions. Therefore, these two layers can be determined first. The other surfaces can also be obtained by limiting the search region based on the previous layers. A 2-D average filter is applied to the segmentation results to obtain smooth surfaces. Finally, the surfaces of ILM, nerve fiber layer (NFL)/ganglion cell layer (GCL), inner plexiform layer (IPL)/inner nuclear layer (INL), INL/outer plexiform layer (OPL), OPL/outer nuclear layer (ONL), outer segment (IS)/outer segment (OS), and RPE/choroid [Fig. 2(c)] are detected by the automatic algorithm.
Fig. 2
Three-dimensional segmentation of the retina by using a 3-D expansion of dynamic programming method. (a) 3-D volume rendering of the retina centered on the macular fovea. (b) Three-dimensional rendering of the seven segmented surfaces from the ILM layer to the RPE/choroid boundary. (c) The corresponding cross-section images overlaid with segmented boundaries indicated by the planes in (a) and (b). The boundaries from inner to outer retina: 1. ILM; 2. NFL/GCL; 3. IPL/INL; 4. INL/OPL; 5. OPL/ONL; 6. IS/OS; and 7. RPE/choroid. Scale bar: 0.4 mm. Video 1 gives a flying view of the segmented layers along the -axis (Video 1, MEPG, 3.18 MB) [URL: http://dx.doi.org/10.1117/1.JBO.20.7.076003.1].

Combined with the 3-D surfaces data, the volumetric microvascular signal of the retina can be segmented into different layers from the IBDV dataset. The en face projection view of the microvascular network for intraretinal layers is produced by using the maximum intensity projection (MIP) method.
3.
Results and Discussion
The 3-D OCT data of the retina from three right eyes was collected to test the feasibility of the method. The seven surfaces of the intraretinal layers were segmented automatically by the developed automated segmentation algorithm. The 3-D volume rendering of the retina centered on the macular fovea with an image area of is shown in Fig. 2(a). The seven segmented surfaces from the ILM layer to the RPE/choroid boundary are reconstructed [Fig. 2(b)] such that the retina from ILM layer to RPE/choroid boundary can be divided into six layers, including NFL, GCL + IPL, INL, OPL, ONL + IS, and OS. Figure 2(c) shows a cross-section image overlaid with the segmented boundaries, and the position of the image is indicated by the planes in Figs. 2(a) and 2(b). The flying view of the segmented layers along the axis is also shown in Video 1. The position of the segmented layers accurately fall on the corresponding boundaries of intraretinal layers [Fig. 2(c) and Video 1], and the six layers of the retina are successfully obtained. Compared with the 2-D image segmentation algorithm in the retina,22,23 the 3-D expansion dynamic programming method can use the information from neighboring B-scans to improve the accuracy and robustness of the algorithm. Garvin et al.24 proposed a graph cut-based method for simultaneous segmentation of multiple 3-D surfaces in the retina. However, the computational complexity was significantly increased in the 3-D data. Compared with the graph cut method, the computational complexity can be significantly reduced by using a dynamic programming method. The processing time for finding an optimal surface from volume data with voxels was with parallel techniques in a Quad-Core Intel i7 2.66GHz workstation. The total processing time for extracting all of the surfaces was and the most time-consuming step was the 3-D median filter processing.
To evaluate the accuracy of the developed automated segmentation algorithm, the results of 27 slides from three sets of 3-D volume data were compared to those achieved by manual segmentation. Nine linear spaced slides (in the fast scan plane) from each dataset were selected, and the fifth slide was centered at the fovea. The manual boundary detections were obtained by three trained graders using ImageJ software [National Institutes of Health (NIH), Bethesda, Maryland]. The automatic detections were extracted from the corresponding positions of the seven segmented surfaces. The deviation of the difference in retinal boundary positions between the automated and manual detections was calculated.23 The boundaries of retinal layers detected by the automated algorithm were very close to the manual segmentation results in both normal and high-myopia eyes. The mean differences were within 2.18 pixels (Table 1) which are similar to 2-D-based segmentation algorithms.22,23 The differences between automatic segmentation and averaged manual segmentation were similar to the difference between manual graders. The layers were segmented based on the normal characteristics of retinal structure. The algorithm may fail with obvious abnormal structures, such as central serous chorioretinopathy (serous detachment of the neurosensory retina) and retinitis pigmentosa (disappearance of the OS layer).30 However, the 3-D expansion dynamic programming method is suitable for interactive segmentation in retinal disease because of a significant reduction in the computational complexity.
Table 1
Differences (mean±SD in pixels) in retinal boundary positions of 27 slides from three sets of 3-D volume data. Between-segmenter indicates the mean differences of M1 versus M2, M1 versus M3, and M2 versus M3; auto indicates the results from automated segmentation; M1, M2, and M3 indicate the results from three trained graders; Avg indicates the averaged results from graders. Each pixel is 2.83 μm.
Retinal layer | Between-segmenter | Auto versus M1 | Auto versus M2 | Auto versus M3 | Auto versus Avg |
---|---|---|---|---|---|
ILM | |||||
NFL/GCL | |||||
IPL/INL | |||||
INL/OPL | |||||
OPL/ONL | |||||
IS/OS | |||||
RPE/choroid |
The microvascular network of the macular retina was also obtained by using the IBDV method. Figure 3(a) shows en face retina vessels from the ILM layer to the RPE/choroid boundary by using the MIP method. The capillary network can clearly be visualized around the foveal avascular zone. Figures 3(b) and 3(c) are cross-sectional views of the vessels merged with the intensity images corresponding to the blue and green lines in Fig. 3(a). Video 2 shows the changes of the vessels along the y axis. Most of the vessel signals are located in the inner retina and choroid. There are also some Doppler variance signals in the outer retina which should be avascular. This should relate to the shadow artifact from the upper layers.10,31 Total processing for obtaining the IBDV images from the raw dataset is after acceleration by the GPU (NVIDIA, GeForce GTX 460).
Fig. 3
In vivo microvascular network of macular retina in a healthy subject (). (a) En face maximum intensity projection (MIP) view image of the retinal microvascular network from ILM layer to RPE/choroid boundary; (b) and (c) cross-sectional views of the retina marked by the blue and green lines in (a). The color-encoded vessel images are merged with intensity images (gray scale). Scale bar: 0.4 mm. Video 2 shows flying view of the vessel signal along the -axis (Video 2, MEPG, 8.47 MB) [URL: http://dx.doi.org/10.1117/1.JBO.20.7.076003.2].

Combined with seven surfaces of the intraretinal data, the microvascular network of the macular retina can be divided into six layers. Figures 4(a)–4(f) present an en face MIP microvascular view of the six layers, including NFL, GCL + IPL, INL, OPL, ONL + IS, and OS + RPE. Most of the vessels in the macular region are located in the inner retina, which agree with the known anatomy.32 In Figs. 4(h) and 4(i), some deeper microvessels may be sheltered from upper vessels in the en face view of the total retina. Segmentation of the individual layers can help to separate these vessels. Microvascular images of individual layers may provide some essential information for early diagnosis of some retinal diseases. The shadow artifacts are also observed in the MIP view of the individual layers [Fig. 4(f)]. However, most of the shadow artifacts are generated from the relatively large vessels [Figs. 3(b) and 3(c)]. Compared with simply projecting the maximum values with a range of depth above the RPE/choroid boundary [Fig. 4(i)], the MIP view from the ILM layer to the RPE/choroid boundary [Fig. 4(h)] can help to improve the quality of an OCT angiography image. Because the retinal vessels are concentrated in the inner retinal layers (GCL, IPL, INL, OPL), the contrast of the microvascular images can be improved by projecting the maximum values from the NFL/GCL boundary to the OPL/ONL boundary [Fig. 4(g)] rather than the total retina [Fig. 4(h)].
Fig. 4
Intensity-based Doppler variance method combined with 3-D segmentation of the retina. (a)–(f): En face MIP view microvascular images of the six individual layers: (a) NFL, (b) GCL + IPL, (c) INL, (d) OPL, (e) ONL + IS, and (f) OS + RPE. (g) MIP view microvascular image of the inner retina (from NFL/GCL boundary to OPL/ONL boundary). (h) MIP view microvascular image of the total retina (from ILM layer to RPE/choroid boundary). (i) MIP view microvascular image of the retina with a range of 400-μm depth above the RPE/choroid boundary. Scale bar: 0.5 mm.

4.
Conclusion
In summary, we have demonstrated a 3-D segmentation method with IBDV based on SS-OCT with a center wavelength of 1050 nm. The feasibility of the method was tested on 3-D data of a retina centered on the macular fovea with an image area of . The seven surfaces of the intraretinal layers were successfully segmented automatically by using a 3-D expansion of a dynamic programming method. The IBDV method was used to acquire the microvascular network of the macular retina. The en face MIP microvascular network images of six layers were obtained after combining the 3-D surfaces of the segmented layers. The morphology of the microvascular network for the individual intraretinal layers can be visualized, and the segmentation method can also be used to enhance the contrast of the vascular images. This method has potential for earlier diagnosis and precise monitoring in retinal vascular diseases.
Acknowledgments
This work was supported by the National Institutes of Health under Grant Nos. R01EY-021529, R01HL-125084, R01HL-105215, and P41-EB015890; and the National Major Equipment Program from the Ministry of Science and Technology in Beijing, China, under Grant No. 2012YQ12008004. Dr. Zhongping Chen has a financial interest in OCT Medical Inc., which, however, did not support this work.
References
Biography
Zhongping Chen is a professor of biomedical engineering and the director of the F-OCT Laboratory at the University of California, Irvine. He is a co-founder and the board chairman of OCT Medical Imaging Inc. His research group has pioneered the development of Doppler optical coherence tomography. He has published more than 200 peer-reviewed papers. He is a fellow of the American Institute of Medical and Biological Engineering, a fellow of SPIE, and OSA.