Automated boundary detection of the optic disc and layer segmentation of the peripapillary retina in volumetric structural and angiographic optical coherence tomography

: To improve optic disc boundary detection and peripapillary retinal layer segmentation, we propose an automated approach for structural and angiographic optical coherence tomography. The algorithm was performed on radial cross-sectional B-scans. The disc boundary was detected by searching for the position of Bruch’s membrane opening, and retinal layer boundaries were detected using a dynamic programming-based graph search algorithm on each B-scan without the disc region. A comparison of the disc boundary using our method with that determined by manual delineation showed good accuracy, with an average Dice similarity coefficient  0.90 in healthy eyes and eyes with diabetic retinopathy and glaucoma. The layer segmentation accuracy in the same cases was on average less than one pixel (3.13


Introduction
Optical coherence tomography angiography (OCTA) is a relatively new method for noninvasive, volumetric imaging of the eye. It uses the blood flow-induced signal variation as the intrinsic contrast mechanism to differentiate vasculature from static tissues [1][2][3][4]. This signal variation can be quantified as the decorrelation between repeated cross-sectional Bscans at each position. The en face projection of the maximum decorrelation based on anatomic tissue slabs produces angiograms which are similar to traditional angiography [5]. We have presented en face angiograms that show the individual layers of the retinal capillary plexuses [6]. The ability of this technique to observe individual capillary networks may enable earlier detection of microvasculopathies. Demonstrating these capillary plexuses is dependent on accurate retinal layer segmentation.
Several methods have been proposed to identify and segment the retinal layers using structural optical coherence tomography (OCT) information [7][8][9][10][11]. We previously proposed a method based on directional graph search, which showed accurate segmentation results even in eyes with macular degeneration and diabetic retinopathy (DR) in the macular region [7,8,10].
However, the region around the optic disc presents a special challenge (Fig. 1). The anatomic disruption caused by the optic disc makes segmentation of the peripapillary retinal layers difficult with the methods shown to be efficient and successful for the macular region. In addition, this region presents significant, inter-individual anatomic variation. Clinically, correct identification of the optic disc boundary and peripapillary retinal layer boundaries is critical for accurate quantification of the optics disc perfusion and distinct capillary plexuses, which is especially important for evaluating diseases like glaucoma [12][13][14]. A different approach is thus necessary to identify the disc boundary and then accurately segment its adjacent peripapillary retinal layers. In this study, we propose an automated disc boundary detection and layer segmentation method for both structural OCT and OCTA of the optic disc. Based on radial B-scans centered on the optic disc, this algorithm first identifies the Bruch's membrane opening (BMO) as the optic disc boundary [15]. Then, each radial B-scan (with the disc portion removed) is segmented along laminar retinal structures using a dynamic programming (DP)-based graph search [16]. These automatically delineated disc boundaries are compared with manual tracing in healthy eyes and eyes with DR and glaucoma.

Data acquisition
The OCTA scans were acquired using a commercial spectral domain OCT system (RTVue-XR, Optovue, Inc., Fremont, CA). The system has a center wavelength of 840 nm with a fullwidth half-maximum bandwidth of 45 nm and an axial scan rate of 70 kHz. The scans were 4.5 × 4.5 mm and 1.6 mm in depth (304 × 304 × 512 pixels) centered on the optic disc. Two repeated B-scans were captured at a fixed position before proceeding to the next location. All 608 B-scans in each data cube were acquired in ~2.9 seconds. The split-spectrum amplitudedecorrelation angiography (SSADA) algorithm [3,17] was used to detect blood flow. The two repeated B-scans were averaged to generate the OCT structural image. For each volumetric data set, two volumetric raster scans, one x-fast scan and one y-fast scan, were registered and merged through an orthogonal registration algorithm to reduce motion artifacts [18]. Figure 1 illustrates the challenges of applying our previous graph search segmentation approach to the optic nerve region. The following layers or boundaries are anatomically important: the inner limiting membrane (ILM), nerve fiber layer (NFL), ganglion cell layer (GCL), inner plexiform layer (IPL), inner nuclear layer (INL), outer plexiform layer (OPL), outer nuclear layer (ONL), retinal pigment epithelium (RPE), and Bruch's membrane (BM). The algorithm segmented the boundaries: vitreous/ILM (red), NFL/GCL (green), IPL/INL (blue), INL/IPL (magenta), OPL/ONL (yellow), IS/OS (red), and RPE/BM (green). However, the presence of the optic disc introduced errors in the segmentation, particularly near the optic disc. Because the algorithm works well for areas not traversing the optic disc, we elected to first identify the disc boundary and remove the disc region on B-scans prior to segmentation. To accomplish this, we first generated radial B-scans from the center of the optic disc, which creates B-scans with a consistent relationship between the peripapillary retina and the optic disc. A flowchart of this process is presented in Fig. 2.  The en face projection of the structural OCT in polar coordinates. The disc center was manually selected (green point in the center). Three hundred and sixty B-scans were radially resampled from 1° to 360° (red lines). Each B-scan ended at the edge of the scan (blue box). (B) The B-scan when the angle was 1°. (C) The Bscan when the angle was 45°. The green lines in panels B and C mark the boundary of the disc, and the blue lines represent the distance from the boundary to the disc center.

Radial B-scan generation
Each radial B-scan consisted of A-lines from the center of the optic disc to the boundary of the volumetric structure. The center of the disc was identified manually before radial resampling ( Fig. 3(A)). Three hundred and sixty radial B-scans were generated from 1° to 360°. We defined the first radial B-scan as oriented towards the fovea. The number of A-lines in each radial B-scan varied depending on the distance from the center of the disc to the edge of the scan (Fig. 3). To generate each radial B-scan, the X-Y raster coordinate was converted to radial coordinate based on the calculated angle and distance, by applying 2-dimensional interpolation. Each radial B-scan was then mirror imaged such that the resulting B-scans presented the peripapillary retina starting at the edge of the scan on the left and the center of the optic disc on the right (Figs. 3(B) and 3(C)). Finally, a 3 × 3 pixel median filter was used on each radial B-scan to reduce speckle noise.

Disc boundary detection
The disc boundary was defined as the position of the BMO, which is at the termination of RPE/BM in each radial B-scan. A grader seeded the BMO in the first B-scan (1°). Before segmentation, all B-scans were truncated to the same length (Fig. 4). Next, the lower boundary of the RPE layer (RPE/BM) was segmented based on the seeded end point using DP-based graph search [16] (red line in Fig. 4(B)). In the graph search method, an intensity gradient map ( , ) o G x z in depth along the A-line was first calculated as follows, In x z is the intensity of the pixel and ( , 1) In x z  is the intensity of the previous pixel located along the A-line. The gradient map was normalized to [0,1] according to the maximum and minimum value, where ( , ) P x z is the cost of the shortest path from the first column of the identified rectangular region to the coordinate ( , ) xz , ( , ) G x z is the pixel value in the corresponding gradient map, () di is one of the six directions, and () wi is the weight of each direction. Here, () wi was set based on experience and testing. The shortest path in the search region manifests as small gradient values across the image (visually, a dark line), and in this case, aligns with the RPE/BM. After the RPE/BM of the first radial B-scan was segmented, the RPE/BM of the other 359 radial B-scans were then segmented based on segmentation from the previous B-scan (Figs. 4(B) and 4(C)). In segmenting each radial B-scan, the search area was limited to 10 pixels above and below the previous RPE/BM boundary. After the RPE/BM was segmented, the difference of pixel values in the gradient map before and after any point in the boundary was calculated as follows, 20 20 ( ) ( , ( )) ( , ( )) 21, 22,..., -20 where () Bj is the position of the RPE/BM in A-line j , ( , ( )) G j B j is the value of the corresponding point in the gradient map, and L is the number of A-lines in the segmented RPE/BM boundary. The point with the largest difference Diff was designated as the position of the BMO (Fig. 4(D)). Finally, 360 BMO positions were interpolated to create a smooth disc boundary.

Automated layer segmentation
After identifying the disc boundary on the radial B-scans, we applied DP-based graph search to automatically segment 6 additional boundaries in the peripapillary retina (Fig. 5). In addition to ( , ) G x z , we generated an inverse gradient map . The main boundaries were segmented one by one according to the shortest path which was calculated using Eq. (3), as previously described [10]. After segmentation was completed on all resampled radial B-scans, 2-dimensional interpolation was used to fill any gaps within a boundary.

Comparison and evaluation
Using the method outlined in the previous section, 10 healthy participants, 10 participants with proliferative diabetic retinopathy (PDR), 10 participants with nonproliferative diabetic retinopathy (NPDR), and 10 participants with glaucoma were imaged. An expert clinical examination determined the diagnoses of the eyes. The participants were enrolled after an informed consent in accordance with an Institutional Review Board approved protocol at Oregon Health and Science University. The study was conducted in compliance with the Declaration of Helsinki.
It took on average 38 seconds to process each volume to generate radial B-scans and identify the disc boundary. The automated layer segmentation in each radial B-scan took on average 12 milliseconds. The processing was done on a workstation with an Inter(R) Xeon(R) CPU E3-1226 v3 @ 3.30GHz and 16.0 GB RAM running MATLAB 2014b (Mathworks, Natick, MA).
Based on the segmentation of the structural OCT images, OCT angiograms were segmented into five slabs: vitreous, radial peripapillary capillary plexus (RCRP), superficial vascular complex (SVC), intermediate capillary plexus (ICP), and deep capillary plexus (DCP). The RPCP was defined as the NFL slab. The SVC was defined as the inner 80% of ganglion cell complex (GCC) which included all structures between the ILM and IPL/INL border. It is dominated by the RPCP in the peripapillary region. The ICP was defined as the outer 20% of GCC and the inner 50% of INL. The DCP was defined as the remaining slab internal to the outer boundary of OPL [14,19]. The projection-resolved (PR) OCTA algorithm was applied to all OCTA scans to remove the interference of flow projection artifacts in the deeper plexuses [20]. En face OCT angiograms were generated by projecting the maximum decorrelation within a slab defined by two relevant boundaries. En face structural OCT was generated by projecting the average reflectance value of each A-line.  Figure 6 shows typical segmentation results and en face angiograms from a healthy eye. The disc boundary was accurately delineated based on the position of BMO (Fig. 6(A)), and the seven boundaries were, likewise, accurately segmented (Figs. 6(B) and 6(C)). The vitreous angiogram was the signal above the ILM and showed slight noise (Fig. 6(D)). The RPCP, SVC, ICP, and DCP were generated based on the segmented boundaries and showed representative vascular patterns (Figs. 6(E)-6(H)).
We also tested our method on OCTA of eyes with DR and glaucoma, which frequently have anatomic alteration in the peripapillary retina and make the segmentation challenging. Figure 7 shows an example of an eye with NPDR. Despite RPE and choroidal atropy and hyperreflectance in the sclera layer (Figs. 7(A)-7(C)), the algorithm was able to identify the BMO. Angiograms of the peripapillary retina demonstrated incongruent areas of nonperfusion in the three plexuses (arrows in Figs. 7(E)-7(G)), consistent with our finding in macular scans in DR [6]. In a case of PDR, retinal neovascularization could be clearly identified as flow signal in the vitreous slab ( Fig. 8(D)).
In glaucoma, NFL thinning can result in segmentation errors in that layer. With this algorithm, the boundaries were accurately segmented even in cases of severe NFL loss (Fig.  9). Figure 9(D) demonstrates isolated loss of RPCP in the NFL with preservation of normal capillary perfusion in the layers below.   To evaluate the accuracy of the disc boundary detection objectively, the Dice similarity coefficient (DSC) between disc regions identified by our method and manual delineation was calculated as follows, where a S and m S are the disc areas detected by our method and an expert grader and am S  is the overlapping area between a S and m S . The mean ± standard deviation DSC value of healthy eyes, eyes with DR, and eyes with glaucoma are shown in Table 1. The automated retinal layer segmentation was also compared with the manually segmented reference boundary. Five eyes from each group (healthy, DR and glaucoma) were randomly chosen. From each scan, 20 B-scans were randomly selected for evaluation. Three graders manually delineated five retinal layer boundaries using Center for Ophthalmic Optics & Lasers-Angiography Reading Toolkit (COOL-ART) [10]. The manually segmented boundaries were averaged among the three graders and taken as the gold standard. The segmentation accuracy of our method, as quantified based on the absolute errors to the gold standard (mean ± standard deviation in units of pixels), is shown in Table 2. Compared to 3D graph search [21], our approach showed less errors in healthy and glaucoma eyes ( Table 3). The segmentation accuracy of the two methods was comparable in DR eyes, with relatively large errors (2-3 pixels) in DR eyes that have more severe distortion of retinal layers. Compared to our previous 2D directional graph search [10], the method proposed in this study was more accurate in all groups of eyes, except the ILM boundary in healthy eyes (Table 4).
Of note, the standard deviation of most boundaries was larger than the mean value. This is due to the thickness of the NFL and the large vessels, which made identification of the layers difficult (Fig. 10).

Discussion
Accurately segmented en face angiograms make visualization of the 3-dimensional OCTA information possible. However, segmentation of the peripapillary area is complicated by the unique anatomy of this region. On B-scans where the disc takes up to a third of the image, connecting the segmentation across the disc can lead to inaccuracies. To improve segmentation in the peripapillary region, we developed an approach which identifies the boundaries of the optic disc and removes the optic disc from consideration when performing layer segmentation. Our method performed well in both healthy and pathologic eyes in comparison to manual segmentation. The segmentation accuracy is superior to currently available 3D graph search and our previous 2D directional graph search methods.
While most of the steps in our approach were automated, there were two minor manual operations. These were identification of the disc center on the en face projection of the structural OCT and identification of the position of BMO in the first radial B-scan. In most cases, the center of the optic disc was near the center of the scan. However, in cases where the scan was significantly off center, we could not use the center of the scan as the center of the disc. A method to identify the disc and then calculate the center on the en face structural OCT would eliminate the need for this manual step [11]. Manually identifying the BMO in the first radial B-scan, while not necessary, allowed for the BM segmentation to be more localized and improved the reliability of the approach.
Beyond OCTA, accurate segmentation is critical in determining the thickness of retinal layers. Thickness maps such as those of the NFL and ganglion cell complex have been shown to be informative in the diagnosis and management of diseases like glaucoma [22,23]. Because our segmentation approach is based on radial B-scans, scans which cover a larger area may require smaller increments during the resampling process to prevent undersampling at the edges of the scan.

Conclusion
We have proposed an optic disc boundary detection and retinal layer segmentation method for structural OCT and OCTA of the peripapillary region. The optic disc boundary was detected based on the position of BMO in radially resampled B-scans. Each radial B-scan with the disc region removed was then segmented using a DP-based graph search algorithm. Our method agreed well with manual grading in both healthy and diseased eyes.