Automated segmentation of retinal layer boundaries and capillary plexuses in wide-field optical coherence tomographic angiography

Advances in the retinal layer segmentation of structural optical coherence tomography (OCT) images have allowed the separation of capillary plexuses in OCT angiography (OCTA). With the increased scanning speeds of OCT devices and wider field images (≥10 mm on fast-axis), greater retinal curvature and anatomic variations have introduced new challenges. In this study, we developed a novel automated method to segment seven retinal layer boundaries and two retinal plexuses in wide-field OCTA images. The algorithm was initialized by a series of points forming a guidance point array that estimates the location of retinal layer boundaries. A guided bidirectional graph search method consisting of an improvement of our previous segmentation algorithm was used to search for the precise boundaries. We validated the method on normal and diseased eyes, demonstrating subpixel accuracy for all groups. By allowing independent visualization of the superficial and deep plexuses, this method shows potential for the detection of plexus-specific peripheral vascular abnormalities.


Introduction
Optical coherence tomography (OCT) [1] is an interferometric imaging technology capable of acquiring high resolution, three-dimensional (3D) images of biological tissue such as the retina through non-invasive and non-contact laser scanning. It has been widely used in the diagnosis of ophthalmic diseases, such as glaucoma [2], diabetic retinopathy (DR) [3], and age-related macular degeneration (AMD) [4], by quantifying the thicknesses of relevant slabs. OCT angiography (OCTA) is a novel clinical tool for the early diagnosis of the diseases affecting retinal circulation and assessment of progression. Based on the variation of OCT signals between B-scans at the same position, OCTA can provide depth-resolved flow signals for the microvasculature. Prior studies have proved that slab-based OCTA can improve the visualization and interpretation of OCTA volumes [5][6][7], and a recent study also showed that vascular abnormalities are better visualized by separating the retinal circulation into three vascular layers [8,9]. Therefore, automated segmentation of the retinal layer boundaries is essential to accurately assess anatomic thickness and capillary plexuses. . However, the wide-field OCT poses new challenges to the existing segmentation algorithms. First, SS-OCT systems, using 1050-nm center wavelength lasers, have decreased the axial resolution and back-scattered reflectance contrast compared to those of the spectral domain commercial devices that use 840-nm center wavelength, which reduces the pixels contained within retinal layers as well as the number of features that can be extracted for machine learning segmentation alternatives. Second, due to the large retinal curvature-associated aberration in the wider field of view, the focusing of wide-field OCT is compromised in the peripheral regions. Third, retinal curvature and anatomic variations are increased as the field of view increases. These characteristics make single source path search algorithms (e.g., graph search) prone to local errors that can be propagated further by the search routine.
Previously, we developed a successful segmentation algorithm based on directional graph search for 3 × 3-and 6 × 6-mm scans of the retina [17]. To address the new challenges associated with wide-field scans, we propose here the Guided Bidirectional Graph Search (GB-GS) method, in which an array of points is used to guide the graph search algorithm in two directions to identify the seven retinal boundaries. The method consists of three steps. First, a guidance point array (GPA) was found to represent the approximate positions of the boundaries. Then, a bidirectional graph search was applied on each point contained in the GPA but not included in any previous paths. The shortest paths between each of the other points in the GPA were used to generate the final boundaries.

Data acquisition
The study was approved by an Institutional Review Board/Ethics Committee of Oregon Health & Science University, and informed consent was collected from all participants, in compliance with the Declaration of Helsinki. Volumetric scans of both eyes were acquired by a prototype 200-kHz SS-OCT system with a 1050-nm central wavelength covering the 10mm (fast-axis) × 8-mm (slow-axis) retinal regions. Two repeated B-scans were taken at each of 400 raster positions, and each B-scan was comprised of 850 A-lines. B-scans at the same position were averaged to improve signal-to-noise ratio of the structural OCT. The OCTA data was calculated by the split-spectrum amplitude decorrelation angiography (SSADA) algorithm [31].

Preprocessing
First, we normalized the B-scan and then flattened it using the center of mass as reference to prevent errors caused by significant tissue curvature [17]. Then, we generated gradient maps to emphasize the transitions between retinal layers with different reflectivity. We reduced the speckle noise by applying a median filter (kernel size -width × height: 3 × 3) and a mean filter (kernel size -width × height: 7 × 3) that preserved the continuity of retinal layer boundaries primarily along horizontal direction. Because the boundaries being segmented exhibited two different intensity transition modes ( Fig. 1) (light-to-dark and dark-to-light) [17], we generated two gradient maps, A G representing dark-to-light transitions and B G representing light-to-dark transitions (Eq. (1)) where ( ) , I x z was the OCT reflectance value at position ( ) , x z , M was the length of Ascans in pixels, and N was the width of B-scans in pixels. From the gradient map A G , we retrieved the boundaries between the vitreous and the inner limiting membrane (ILM), the inner nuclear layer (INL) and the outer plexiform layer (OPL), as well as the upper boundary of the ellipsoid zone (EZ) (Fig. 2(A)). From the gradient map A G , we retrieved the remaining four boundaries that were between the nerve fiber layer (NFL) and the ganglion cell layer (GCL), between the inner plexiform layer (IPL) and the inner nuclear layer (INL), between the OPL and the outer nuclear layer (ONL), and between the retinal pigment epithelium (RPE) and Bruch's membrane (BM) (Fig. 2(B)).

Guidance point array
In this step, we generated for each boundary an array of points indicating its approximate position based on information extracted from the gradient maps. This GPA regulates the subsequent bidirectional graph search for the actual layer boundaries. GPAs were generated in a pre-determined order, taking advantage of the characteristics of gradient maps and retinal anatomy to minimize deviations from the correct boundaries (Fig. 3). First, the vitreous/ILM and upper EZ boundaries were processed from gradient map A G , as they exhibited the greatest contrast with surrounding tissue. Then, using the EZ layer as the upper boundary, the set of points corresponding to the RPE/BM's GPA was recognized from the B G gradient map. Subsequently, the upper boundary was fixed at the vitreous/ILM boundary, and B G was used to sequentially extract the GPA for the OPL/ONL, which had the EZ layer as the lower boundary. Then we extracted the GPAs for the IPL/INL and NFL/GCL, for which each GPA served as the lower boundary of the next GPA. Finally, the GPA for the INL/OPL was generated from the gradient map A G using the IPL/INL and OPL/ONL as upper and lower boundaries respectively. The first GPAs to be identified were the vitreous/ILM and upper EZ boundaries, which were not limited by any reference boundaries. To localize them, we first reduced speckle noise by down-sampling the gradient map A G and the B-scan by a factor of five to a size of 170 × 208 pixels. The vitreous/ILM boundary plays a very import role in subsequent operations. For the GPA identification, we compounded a new B-scan with enhanced contrast between the vitreous and the ILM by adding the gradient map A G to the normalized B-scan. We then binarized the enhanced B-scan by thresholding pixels with OCT reflectance values below the average reflectance value, which removes nonzero pixels in the vitreous. Then, the first nonzero pixels in each A-line were selected to form the GPA of the vitreous/ILM.
The second GPA to be recognized was the upper EZ boundary. Either this or the previously identified vitreous/ILM boundary contain the lowest gradient values in each A-line in the map A G , making it easy to identify the EZ. Then, the binary image was up-sampled to the original number of pixels, and the 170 GPA points identified were reassigned to the Alines with indices 5n + 1 (n = 0…169).
After the first two GPAs were generated from enhanced B-scans, the remaining five were obtained from the corresponding gradient maps, searching one of every five A-lines restricted to the corresponding upper and lower boundaries assigned above. These GPAs were first enhanced by a horizontal gradient operator (Eq. (2)), and the first point with parameter t<-0.02 (where t is the threshold assigned to ' B G ) was selected for the corresponding GPA (Fig.  4). Due to the relatively low contrast of image, points contained in the GPA were occasionally distant from the actual boundary (Fig. 5). Based on the prevalence of noise and relative flat GPA curves in the wide field of view, we used a mean filter (kernel size -width × height: 9 × 1) on the GPA to remove unreliable points and ensure the accuracy of the operation described below in Section 2.4.

Guided bidirectional graph search
Once the GPAs were identified, we implemented a guided bidirectional graph search algorithm for retinal layer segmentation (Fig. 6(A)). For any point S, we searched for graph points in two directions (left and right). For the next point L (or R), we appointed 5 nearby candidate points in the adjacent A-line (Figs. 6(B-C)) and chose the one with minimum gradient as the next node in the path. Unlike our previous directional graph search algorithm [17], we started from a virtual point located outside the image (Fig. 6(B)) and crossed a collection of points that may or may not fall in the GPA of the boundary under scrutiny. All GPA points crossed by this searched path were dropped from subsequent analysis (Figs. 7(A-B)), and guided bidirectional graph search started again from the next GPA point that was not contained in any previous graph recognized for the current boundary ( Fig. 7(B)). This process was repeated, generating each time a potentially different graph until all GPA points belonged to one of the graphs (Figs. 7(B-E)). Then, all graphs thus generated were merged by the rationale explained in Section 2.5 below.  . GPA points that were located on the first path (red points) and points left out of the graph (blue asterisks) were identified, and a second path was created bi-directionally by graph search, starting from the first GPA point left out of the previous path (in B, blue star, red arrow). Blue asterisks crossed by the second path became red points and did not trigger the start of a future graph search. The process was repeated (C-E) until all blue asterisks eventually form part of one candidate path.
Although there were enough points in the GPA to support a point-to-point shortest path search, we preferred the bidirectional graph search to detect the boundary because we observed that some points in the GPA were outside the manually-segmented interfaces. Therefore, the graph of the layer boundary should not be forced to cross all GPA points, and a different boundary detection and merging scheme was necessary.

Path merging
The preceding procedures generated several possible paths for each boundary in a B-scan. To obtain the final boundary, we evaluated the deviation of each candidate path from the GPA in sections of a B-scan (Fig. 8). For example, from an interval bounded by three points of the GPA with indices a, b, and c, we selected the most accurate of all paths within this interval and assigned it to all A-lines with indices between a and b. To decide the most accurate path within an interval, we designed the evaluation function in (Eq. (3)).  8). Then, the process was repeated for the A-lines in the following interval, i.e., with indices between b and c, etc.  (3)). Two intervals between GPA points a, b, and c were emphasized, and three different paths similar to those generated in (Fig. 7) were represented in light blue, dark blue, and orange color. According to (Eq. (3)), the pixels crossed by the light blue path were assigned to the final path between points a and b, and the pixels crossed by the dark blue path were assigned to the final path between b and c.

Segmentation of capillary plexuses
We extracted two vascular plexuses from the segmented OCTA volume (Figs. 9(A-B)): the superficial vascular complex (SVC) (Fig. 9(C) and the deep vascular complex (DVC) (Fig.  9(D)). En face angiograms of the capillary plexuses were generated by the maximum projection of OCTA flow signals within the slab.

Study population
We tested our segmentation method on normal eyes and eyes with glaucoma, diabetic retinopathy, and retinitis pigmentosa (Table 1). For all cases the seven layers were segmented to identify the vitreous/ILM, NFL/GCL, IPL/INL, INL/OPL, OPL/ONL, EZ, and RPE/BM.

Segmentation performance
We ran the GB-GS algorithm in Matlab R2017a on a desktop PC equipped with an Intel(R) Core(TM) i7-6700K @4.0GHz CPU and 32 GB RAM. The average run time of our algorithm was 0.3 seconds per B-scan. Our method correctly segmented retinal layer boundaries, even in the areas of large vessel shadows (Fig. 10(A)) and small cysts (Figs. 10(B-C)). Segmentation errors were present in areas of extremely low contrast between layers ( Fig.  10(D)); in areas with retinal neovascularization, which could significantly affect the surface of the ILM (Fig. 10(E)); and in an area with a partially separated epiretinal membrane ( Fig.  10(F)). To evaluate segmentation accuracy, we compared the automatic segmentation results with manual segmentation performed by a masked grader. For each eye, 20 B-scans of one volumetric data set were randomly selected for evaluation. The position of the manual boundary was subtracted from the position of the automatic boundary without any manual corrections in all A-lines under scrutiny, and the segmentation accuracy was determined ( Table 2). Subpixel accuracy was present for the four groups, with the most accurate being the vitreous/ILM boundary, which is the one with the highest perceived contrast. 5.12 ± 13.64 Differences in segmentation between manual grading and automated grading were measured in micron and presented as means ± standard deviations. The digital pixel size in the axial direction was 5.5 µm.
Thanks to the stability and robustness of GB-GS, our method can also be used to segment small field of view OCT scans (3 × 3-and 6 × 6-mm). To evaluate the performance on these images, we randomly selected 20 volumetric scans acquired by a 70-kHz commercial AngioVue system (RTVue-XR; Optovue, Inc.) (  We randomly selected 31 B-scans from each volumetric scan, for a total of 620 B-scans. For each B-scan, we applied these three methods to segment 7 retinal boundaries, respectively. The segmentation results of the three methods were compared to manual grading (Table 4). Apparently, our method is superior to other two methods on at least five of seven layers.

Clinical applications
To evaluate the benefits of our segmentation method in the computation of clinically useful parameters, we applied it to the detection of the non-perfusion area in one eye with DR. Capillary nonperfusion is an important feature of DR [34,35], and quantification of it may be an important biomarker of disease progression. In particular, the larger scanning area of widefield OCTA will likely improve the sensitivity of this metric for early stages of the disease because the manifestations of capillary dropout in DR begin in the peripheral retina rather than the central macula.
Using our automated segmentation method, we segmented each layer on a structural OCT B-scan (Fig. 11(A)). The en face angiogram of the SVC and DVC flow were generated (Figs. 11(B-C)), and a slab subtraction algorithm [6,7,36] was applied to reduce the prevalence of projection artifacts in the DVC. Then, we generated a nonperfusion map (Fig. 11(D)) using an automated algorithm developed previously [6,7]. The resulting images demonstrated areas of capillary nonperfusion over 7.04 mm 2 that were specific to individual plexuses (Figs. 11(B-C)), allowing plexus-specific detection of nonperfusion in OCTA. Another possible use of wide-field OCTA is identification of neovascularization in DR eyes. A 10 × 25-mm wide-field OCTA, produced by montaging four scans, demonstrates a large area of neovascularization temporal to the macula (Fig. 12). Because wide-field-OCTA visualizes the neovascularization clearly without leakage, quantification of neovascularization is possible, allowing objective monitoring of treatment response. Fig. 12. Wide-field OCTA of a patient with proliferative diabetic retinopathy. A large area of neovascularization (yellow) temporal to the macula was present. This image is montaged from four 10 × 8-mm scans. The total size is 10 × 25-mm. The traditional 3 × 3-and 6 × 6-mm commercial OCTA images at the central macular area are indicated by dashed squares respectively. Unlike the fluorescein angiograms, OCTA demonstrates the neovascularization clearly without leakage and allows for quantification.

Discussion
In this study, we demonstrated an improvement over our previous graph search retinal layer segmentation algorithm and OCTExplorer algorithm to achieve a more accurate delineation of the seven layer boundaries imaged by wide-field OCT scanning. The method was able to segment both healthy and diseased retinas, including hypo-reflective regions affected by vascular shadows and retinal cysts.
The main advantage of this algorithm is the ability to accurately segment retinal layers over a large scanning area. Traditional OCTA had been restricted from its inception to narrow fields of view, i.e., 3 × 3-, 4.5 × 4.5-, and 6 × 6-mm, which are still standard in commercial machines. Wide-field OCTA is a natural evolution of this technology, compelled by the clinical demand for better visualization of the peripheral retina. Stitching many images by registration techniques is an alternative to generate retinal angiograms of larger size, and it is inherently better to montage a few wide-field scans (e.g., 10 × 6-mm) than numerous narrowfield scans (e.g., 6 × 6-mm). For instance, the angiogram represented in Fig. 12 was generated by montaging of four 10 × 8-mm scans, whereas at least ten 6 × 6-mm scans would be needed to represent the same area. However, the advantage of wide-field scanning comes at the expense of more challenging segmentation across ultra-wide B-scans. Our method based on GB-GS not only can handle the macular area, but also can accurately segment the optic disc region and peripheral retinal region.
Recently, segmentation of retinal layers and pathological structures has also been accomplished by alternative supervised machine learning methods such as deep learning [21,22]. An advantage of our current guided graph search method is that unlike deep learning solutions, it does not need a large, annotated data set to be used for network training, and hence it is suitable for small studies, for data acquired by lab-built prototype devices, and for diseases in which even manual segmentation of boundaries is uncertain and could introduce confusion during training. Moreover, the machine learning methods reported previously only generated probability maps and still needed a post-processing step (e.g., graph search or conditional random fields) to generate sharp boundaries. In contrast, our results show that the method proposed here is generalizable to different retinal pathologies. This method is superior to previous graph search solutions in that it considers the laminar structure of the retina and performs the search in two directions, relying on the GPA to prevent graph deviations from the anatomically connected boundaries. Finally, segmentation is performed faster than machine learning alternatives owing to the lower computational requirements.
The limitations of the software can be summarized as follows. First, the method depends strongly on the gradient information at the layer boundaries and might fail for acquisitions with extremely low contrast between layers. Second, due to the order in which boundary detection is defined, the segmentation is sensitive to any errors in segmentation of previous graphs bounding the position of its upper and lower limits. To address this issue, the boundaries least likely to be erroneously segmented owing to the highest contrast were chosen to precede the segmentation of the boundaries more likely to be affected by disease.

Conclusions
We proposed a novel automatic segmentation method to find the boundaries of seven retinal layer boundaries in wide-field OCTA images. Our algorithm showed sub-pixel accuracy in both normal and diseased eyes. The extraction of thin slab boundaries over a large area has great potential for use in the improved diagnosis and progression assessment of diseases. This is especially true for diseases that begin from the peripheral retina and affect large areas, such as DR and inherited retinal diseases, where evaluation by OCTA was limited in the past to a small field of view.