Shared-hole graph search with adaptive constraints for 3 D optic nerve head optical coherence tomography image segmentation

Optic nerve head (ONH) is a crucial region for glaucoma detection and tracking based on spectral domain optical coherence tomography (SD-OCT) images. In this region, the existence of a “hole” structure makes retinal layer segmentation and analysis very challenging. To improve retinal layer segmentation, we propose a 3D method for ONH centered SD-OCT image segmentation, which is based on a modified graph search algorithm with a shared-hole and locally adaptive constraints. With the proposed method, both the optic disc boundary and nine retinal surfaces can be accurately segmented in SD-OCT images. An overall mean unsigned border positioning error of 7.27 ± 5.40 μm was achieved for layer segmentation, and a mean Dice coefficient of 0.925 ± 0.03 was achieved for optic disc region detection. © 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement OCIS codes: (100.0100) Image processing; (100.2960) Image analysis; (170.4470) Ophthalmology; (170.4500) Optical coherence tomography. References and links 1. M. R. Hee, J. A. Izatt, E. A. Swanson, D. Huang, J. S. Schuman, C. P. Lin, C. A. Puliafito, and J. G. Fujimoto, “Optical coherence tomography of the human retina,” Arch. Ophthalmol. 113(3), 325–332 (1995). 2. G. J. Jaffe and J. Caprioli, “Optical coherence tomography to detect and manage retinal disease and glaucoma,” Am. J. Ophthalmol. 137(1), 156–169 (2004). 3. Y. C. Tham, X. Li, T. Y. Wong, H. A. Quigley, T. Aung, and C. Y. Cheng, “Global prevalence of glaucoma and projections of glaucoma burden through 2040: a systematic review and meta-analysis,” Ophthalmology 121(11), 2081–2090 (2014). 4. R. N. Weinreb, T. Aung, and F. A. Medeiros, “The pathophysiology and treatment of glaucoma: a review,” JAMA 311(18), 1901–1911 (2014). 5. H. A. Quigley and A. T. Broman, “The number of people with glaucoma worldwide in 2010 and 2020,” Br. J. Ophthalmol. 90(3), 262–267 (2006). 6. H. A. Quigley, “Number of people with glaucoma worldwide,” Br. J. Ophthalmol. 80(5), 389–393 (1996). 7. R. Bock, J. Meier, L. G. Nyúl, J. Hornegger, and G. Michelson, “Glaucoma risk index: automated glaucoma detection from color fundus images,” Med. Image Anal. 14(3), 471–481 (2010). 8. J. Nayak, R. Acharya U, P. S. Bhat, N. Shetty, and T. C. Lim, “Automated diagnosis of glaucoma using digital fundus images,” J. Med. Syst. 33(5), 337–346 (2009). 9. H. Ahmad, A. Yamin, A. Shakeel, S. O. Gillani, and U. Ansari, “Detection of glaucoma using retinal fundus images,” in Robotics and Emerging Allied Technologies in Engineering (iCREATE), 2014 International Conference on. (IEEE, 2014), pp. 321–324. 10. R. Bock, J. Meier, G. Michelson, L. G. Nyul, and J. Hornegger, “Classifying glaucoma with image-based features from fundus photographs,” in Joint Pattern Recognition Symposium. (Springer, 2007), pp. 355–364. 11. N. Inoue, K. Yanashima, K. Magatani, and T. A. K. T. Kurihara, “Development of a simple diagnostic method for the glaucoma using ocular Fundus pictures,” in 27th Annual International Conference. Engineering in Medicine and Biology Society, 2005. (IEEE-EMBS, 2005), PP. 3355–3358. 12. P. S. Mittapalli and G. B. Kande, “Segmentation of optic disk and optic cup from digital fundus images for the assessment of glaucoma,” Biomed. Signal Process. Control 24, 34–46 (2016). Vol. 9, No. 3 | 1 Mar 2018 | BIOMEDICAL OPTICS EXPRESS 962 #312385 https://doi.org/10.1364/BOE.9.000962 Journal © 2018 Received 13 Nov 2017; revised 8 Jan 2018; accepted 23 Jan 2018; published 2 Feb 2018 13. L. Chakrabarty, G. D. Joshi, A. Chakravarty, G. V. Raman, S. R. Krishnadas, and J. Sivaswamy, “Automated Detection of Glaucoma From Topographic Features of the Optic Nerve Head in Color Fundus Photographs,” J. Glaucoma 25(7), 590–597 (2016). 14. J. C. Mwanza, J. D. Oakley, D. L. Budenz, and D. R. Anderson; Cirrus Optical Coherence Tomography Normative Database Study Group, “Ability of cirrus HD-OCT optic nerve head parameters to discriminate normal from glaucomatous eyes,” Ophthalmology 118(2), 241–248 (2011). 15. B. C. Chauhan and C. F. Burgoyne, “From clinical examination of the optic disc to clinical assessment of the optic nerve head: a paradigm change,” Am. J. Ophthalmol. 156(2), 218–227 (2013). 16. J. C. Downs, H. Yang, C. Girkin, L. Sakata, A. Bellezza, H. Thompson, and C. F. Burgoyne, “Three-dimensional histomorphometry of the normal and early glaucomatous monkey optic nerve head: neural canal and subarachnoid space architecture,” Invest. Ophthalmol. Vis. Sci. 48(7), 3195–3208 (2007). 17. N. G. Strouthidis, H. Yang, J. C. Downs, and C. F. Burgoyne, “Comparison of clinical and three-dimensional histomorphometric optic disc margin anatomy,” Invest. Ophthalmol. Vis. Sci. 50(5), 2165–2174 (2009). 18. A. S. Reis, N. O’Leary, H. Yang, G. P. Sharpe, M. T. Nicolela, C. F. Burgoyne, and B. C. Chauhan, “Influence of Clinically Invisible, but Optical Coherence Tomography Detected, Optic Disc Margin Anatomy on Neuroretinal Rim Evaluation Clinically Invisible Optic Disc Margin Anatomy,” Invest. Ophthalmol. Vis. Sci. 53(4), 1852–1860 (2012). 19. B. J. Antony, M. S. Miri, M. D. Abràmoff, Y. H. Kwon, and M. K. Garvin, “Automated 3D segmentation of multiple surfaces with a shared hole: segmentation of the neural canal opening in SD-OCT volumes,” in International Conference on Medical Image Computing and Computer-Assisted Intervention. (Springer, 2014), pp. 739–746. 20. S. J. Chiu, X. T. Li, P. Nicholas, C. A. Toth, J. A. Izatt, and S. Farsiu, “Automatic segmentation of seven retinal layers in SDOCT images congruent with expert manual segmentation,” Opt. Express 18(18), 19413–19428 (2010). 21. A. Lang, A. Carass, M. Hauser, E. S. Sotirchos, P. A. Calabresi, H. S. Ying, and J. L. Prince, “Retinal layer segmentation of macular OCT images using boundary classification,” Biomed. Opt. Express 4(7), 1133–1152 (2013). 22. Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” Nature 521(7553), 436–444 (2015). 23. S. Apostolopoulos, S. De Zanet, C. Ciller, S. Wolf, and R. Sznitman, “Pathological OCT Retinal Layer Segmentation using Branch Residual U-shape Networks,” in International Conference on Medical Image Computing and Computer-Assisted Intervention, (September, 2017), pp. 294–301. 24. L. Fang, D. Cunefare, C. Wang, R. H. Guymer, S. Li, and S. Farsiu, “Automatic segmentation of nine retinal layer boundaries in OCT images of non-exudative AMD patients using deep learning and graph search,” Biomed. Opt. Express 8(5), 2732–2744 (2017). 25. Z. Hu, M. Niemeijer, K. Lee, M. D. Abramoff, M. Sonka, and M. K. Garvin, “Automated segmentation of the optic disc margin in 3-D optical coherence tomography images using a graph-theoretic approach,” Proc. SPIE 7262, 72620U (2009). 26. Z. Hu, C. A. Girkin, A. Hariri, and S. R. Sadda, “Three-dimensional choroidal segmentation in spectral OCT volumes using optic disc prior information,” Proc. SPIE 9697, 96971S (2016). 27. B. J. Antony, M. D. Abràmoff, K. Lee, P. Sonkova, P. Gupta, Y. Kwon, M. Niemeijer, Z. Hu, and M. K. Garvin, “Automated 3D segmentation of intraretinal layers from optic nerve head optical coherence tomography images,” Proc. SPIE 7626, 76260U (2010). 28. K. Lee, M. Niemeijer, M. K. Garvin, Y. H. Kwon, M. Sonka, and M. D. Abramoff, “Segmentation of the optic disc in 3-D OCT scans of the optic nerve head,” IEEE Trans. Med. Imaging 29(1), 159–168 (2010). 29. F. Shi, B. Tian, W. Zhu, D. Xiang, L. Zhou, H. Xu, and X. Chen, “Automated choroid segmentation in threedimensional 1-μm wide-view OCT images with gradient and regional costs,” J. Biomed. Opt. 21(12), 126017 (2016). 30. P. Zang, S. S. Gao, T. S. Hwang, C. J. Flaxel, D. J. Wilson, J. C. Morrison, D. Huang, D. Li, and Y. Jia, “Automated boundary detection of the optic disc and layer segmentation of the peripapillary retina in volumetric structural and angiographic optical coherence tomography,” Biomed. Opt. Express 8(3), 1306–1318 (2017). 31. E. Gao, F. Shi, W. Zhu, C. Jin, M. Sun, H. Chen, and X. Chen, “Graph Search–Active Appearance Model based Automated Segmentation of Retinal Layers for Optic Nerve Head Centered OCT Images,” in SPIE Medical Imaging, (SPIE, 2017), pp. 101331Q–101331Q. 32. L. Breiman, “Random forests,” Mach. Learn. 45(1), 5–32 (2001). 33. G. Staurenghi, S. Sadda, U. Chakravarthy, and R. F. Spaide; International Nomenclature for Optical Coherence Tomography (IN•OCT) Panel, “Proposed lexicon for anatomic landmarks in normal posterior segment spectraldomain optical coherence tomography: the IN•OCT consensus,” Ophthalmology 121(8), 1572–1578 (2014). 34. K. Li, X. Wu, D. Z. Chen, and M. Sonka, “Optimal surface segmentation in volumetric images-a graph-theoretic approach,” IEEE Trans. Pattern Anal. Mach. Intell. 28(1), 119–134 (2006). 35. M. K. Garvin, M. D. Abràmoff, X. Wu, S. R. Russell, T. L. Burns, and M. Sonka, “Automated 3-D intraretinal layer segmentation of macular spectral-domain optical coherence tomography images,” IEEE Trans. Med. Imaging 28(9), 1436–1447 (2009). 36. M. K. Garvin, M. D. Abràmoff, R. Kardon, S. R. Russell, X. Wu, and M. Sonka, “Intraretinal layer segmentation of macular optical coherence tomography images using optimal 3-D graph search,” IEEE Trans. Med. Imaging 27(10), 1495–1505 (2008). Vol. 9, No. 3 | 1 Mar 2018 | BIOMEDICAL OPTICS EXPRESS 963


Introduction
Optical coherence tomography (OCT) is a noninvasive imaging technology that produces high-resolution cross-sectional images of retina [1].It has been increasingly used for diagnosing and managing a variety of eye diseases, such as glaucoma, macular hole, agerelated macular degeneration (AMD), and diabetic macular edema (DME) [2].Among them, glaucoma is the leading cause of irreversible blindness globally.The number of people with glaucoma worldwide will increase to 111.8 million in 2040 based on a report [3], disproportionally affecting Asian and African people.Because glaucoma may be asymptomatic until a relatively late stage, diagnosis is often delayed [4].As a result, improved methods for glaucoma screening and quantitative analysis are urgently needed [5,6].
Most existing methods for glaucoma disease analysis were based on fundus images in two dimensions [7][8][9][10][11][12][13]. 2D features were incorporated into their methods, but the 3D structural information is missing.With the introduction of spectral domain optical coherence tomography (SD-OCT), high resolution imaging of the 3D structure of optic nerve head (ONH) becomes possible.The 3D structural information such as retinal layer thickness and optic disc boundary position has the excellent ability to discriminate normal eyes and mild glaucoma eyes.With the increasing of intraocular pressure, the retinal nerve fiber layer (RNFL) becomes thinner.According to [14], the comparison between normal and eyes with even mild glaucoma showed significant differences in the thickness of the RNFL.On the other hand, the optic disc boundary, also called Bruch's membrane opening (BMO) showed in ONH-centered OCT image is a stability landmark, as it is unaffected by glaucomatous cupping [15].The BMO can be used for computing ONH parameters, and tracking the glaucomatous changing [16][17][18].Therefore, a reliable SD-OCT image analysis method for optic disc boundary detection and retinal layer segmentation in the ONH region is urgently needed for glaucoma disease evaluation and treatment planning.
B-scans of ONH centered SD-OCT image obtained from a normal eye and an eye with glaucoma is showed in Fig. 1, with optic disc boundary marked using dotted yellow lines.In both cases, the retinal structure changes dramatically near ONH, and no retinal layers exist inside the optic disc region.Therefore, the retinal layers in ONH centered SD-OCT image have a shared hole.Such discontinuity makes it difficult to accurately segment the layers when the exact position of the hole boundary, i.e., the optic disc boundary are not known a priori [19].For decades, several studies were devoted to OCT image layer segmentation and analysis.Chiu et al. [20] proposed an automatic method for segmenting retinal layers using graph theory and dynamic programming (GTDP).Lang et al. [21] built a random forest classifier to segment eight retinal layers in macular centered OCT images.Recently, deep learning [22] has been demonstrated to be a powerful tool in many fields.For OCT image layer segmentation and analysis, Apostolopoulos et al. [23] proposed a fully-convolutional convolutional neural networks (CNN) architecture for OCT retinal layer segmentation, which combined dilated residual blocks in an asymmetric U-shape configuration.Fang et al. [24] developed a framework combining CNN and graph search methods for nine retinal layers segmentation on OCT images with non-exudative AMD, the CNN based probability maps were used to create the final boundaries by GTDP method.
For ONH centered OCT image layer segmentation and analysis, Hu et al. [25] developed a graph theory based approach for optic disc segmentation in OCT images.First, the conventional graph search method was used to segment the Bruch's membrane.Then, a projection image was generated using the Bruch's membrane as reference.Then the optic disc boundary was detected by 2D graph search method in the projection image.In [26], a similar method was used to detect the optic disc region and the graph search constraints were modified accordingly to achieve automatic choroid segmentation in ONH-centered SD-OCT images.In [27], Antony et al. used graph search method to segment 7 surfaces (6 layers) from ONH centered SD-OCT images with statistical smoothness constraints.The RNFL thickness of normal and glaucomatous retina was compared.Lee et al. [28] segmented four intra-retinal surfaces and used their positions as reference to extract several features for optic cup and rim classification.Shi et al. [29] approximately detected the optic disc boundary in wide-view SS-OCT images, by finding the discontinuity of the inner-outer retina boundary based on modified Canny operator.The result was used for subsequent choroid segmentation.Zang et al. [30] detected the optic disc boundary by searching for the position of Bruch's membrane opening firstly, and then the retinal layer boundaries were detected with a dynamic programming-based graph search algorithm as shortest path in each OCT B-scan.Gao et al. [31] combined the active appearance model (AAM) and graph search method (GS-AAM) for ONH centered OCT images layer segmentation.The global shape constraints learned from the statistical model were used to constrain the graph search method.
All of the methods mentioned above for ONH centered OCT image analysis did not explicitly consider the hole structure when applying graph search methods, although some of them modified the interrelation constraints or the cost function based on the location of optic disc.Ignoring the existence of the optic disc region will make layer segmentation prone to error near the ONH region.This is because the deep structure inside the optic disc region is different than the other regions and will mislead the segmented surfaces away from their actual locations near the boundary of the hole.On the other hand, detection of the optic disc boundary is challenging due to the existence of large vessels and the presence of external oblique border tissues attaching to the end of BM surface.The errors are especially large in the B-scans where the optic disc boundaries are close.In this paper, we propose a shared-hole graph search method with locally adaptive constraints for multiple surfaces segmentation in ONH centered SD-OCT images.After preprocessing, the volumetric image is first polar-transformed, so that the optic disc boundaries are kept apart.Graph search method with locally adaptive constraints is first used to detect four surfaces in the polar-transformed image.Then the random forests classifier [32] is used to detect the optic disc region with several novel surface-derived features.With the position of optic disc known, new graphs are constructed with particular arcs specifying the shared-hole.For accurate and robust layer segmentation, locally adaptive surface smoothness constraints are learned from previous segmented layers.Finally, 9 retinal surfaces and optic disc boundary can be obtained.Figure 2 shows an ONH centered SD-OCT B-scan of normal eye with the 9 surfaces that define 8 retinal layers or layer complexes.The surfaces are numbered 1 to 9 from top to bottom, and the retinal layers are defined as nerve fiber layer (NFL), ganglion cell layer & inner plexiform layer (GCL + IPL), inner nuclear layer (INL), outer plexiform layer (OPL), outer nuclear layer & henle fiber layer (ONL + HFL), ellipsoid zone (EZ), outer segments of photoreceptors (OSP), interdigitation zone & RPE/Bruch's complex (IZ + RPE) [33].
The contributions of this paper are summarized as follows: (1) An accurate optic disc boundary detection method is proposed based on features representing local intensity, structure and texture information, and is effective both for normal and glaucomatous data.
(2) A novel graph construction method suitable for structure with a shared-hole is proposed, which improves segmentation accuracy near the optic disc boundary.
(3) A graph search method with locally adaptive constraints is proposed for multi-surface segmentation.The varying constraints can adapt to smoothness changes of surfaces, and thus both flat and steep surfaces can be segmented accurately in ONH centered OCT images.

Methods
In this section, we first review the conventional graph search method, and then introduce our proposed shared-hole graph search method with locally adaptive constraints for accurate layer segmentation.Then the proposed method is described in details, which consists of three parts: pre-processing, optic disc boundary detection and surface segmentation.The pre-processing includes denoising, B-scan alignment and polar transformation.In the step of optic disc boundary detection, a supervised random forests classifier with novel surface-derived features is applied on the image columns, and then a 2D graph search method is used to refine the optic disc boundary detected by the classifier.Based on the optic disc boundary detection result, a novel shared-hole graph search method with locally adaptive constraints is proposed to segment 9 retinal surfaces in the original B-scans.The flowchart of our proposed method is shown in Fig. 3.

Review of graph search
The graph search algorithm originally proposed by Li et al. [34] was designed to achieve globally optimal surface segmentation in volumetric data such as medical images.It was successfully applied in retinal layer segmentation [35][36][37][38][39].In the graph search algorithm, volumetric image is defined as a 3D matrix ( , , ) I x y z with size X Y Z × × , and the surface is defined by a function ( , ) z f x y = , where For single surface segmentation, two smoothness constraints, x Δ and y Δ are enforced to ;avoid the change of surface height larger than certain range between neighboring surface points along x and y-direction, respectively.The cost function ( , , ) c x y z is allocated for each voxel in the volumetric image.This cost is inverse to the likelihood of each voxel belongs to the desired surface, so that the optimal surface is the one with the minimum cost.A nodeweighted directed graph ( , ) G V E is constructed from the volumetric image.Each node in V corresponds to one and only one voxel in ( , , ) I x y z .The weight of each node is defined as ( , , ) z = 0 ( , , ) ( , , ) ( , , 1) otherwise c x y z w x y z c x y z c x y z The arc set E consists of intra-column arcs and inter-column arcs.The intra-column arcs connect each node with its immediate neighbor below with infinite cost, and the inter-column arcs are constructed according to the smoothness constraints.For multiple surfaces segmentation, the desired multiple surfaces are constrained not only by smoothness constraints of each surface, but also by surface separation constraints specifying their interrelations.K sub-graph ( , ) ( 1,...., ) surfaces segmentation.Inter-surface arcs are enforced to restrict the maximum and minimum horizontal distances between surface pairs.Figure 4 In general, the goal of graph search method is to find the optimal surfaces which satisfy the following conditions.(1) Each surface satisfies its own surface smoothness constraints.
(2) Each pair of surfaces satisfies their surface separation constraints.(3) The total cost of voxels on surfaces is minimized.This cost minimization problem is then transformed into a Min-Cut/Max-Flow problem [40].
The major limitation of the conventional graph search algorithm lies in the fact that the detected optimal surface must be continuous in the whole volumetric data.But for cases such as ONH centered OCT image segmentation, where a hole exists on the surfaces, the conventional graph search method may fail to give satisfying results.The surfaces may be misled by deeper structures inside the ONH, causing segmentation errors at the hole boundary.To solve this problem, the graph construction needs to be modified to enforce appropriate constraints on the boundary points.
Another limitation of conventional graph search is that the surface smoothness constraints are constant for all locations.However, the surfaces may be less smooth in certain locations than the others.For example, the upper surface of NFL drops suddenly near the optic disc region.Using a large global smoothness constraint will make the detected surface positions more easily affected by noise and other artifacts.Therefore, locally adaptive constraints are needed for accurate surface segmentation.
In [26,27,31,35,41] varying constraint graph search was used for retinal layer segmentation.The constraints were computed from statistic models learned from training data.As a supervised method, it worked well for normal retinas, but may get into trouble when applied to retinas with different type of diseases, and to complicated and changeable ONH centered SD-OCT image.To solve the deficiency of conventional graph search method, we propose shared-hole graph search method for single and multiple surfaces segmentation.The key innovation of the proposed shared-hole graph search method is its novel graph construction.Special hard constraint arcs prevent the surface being misled by deep structures in the optic disc region, while maintaining the segmentation problem as a minimum closed set problem solvable with global optimization.

Shared-hole graph search
The graph construction for single surface segmentation is illustrated in Fig. 5.Each red node represents one and only one voxel in ONH centered SD-OCT image.The translucent yellow region indicates the detected optic disc region.The blue arcs represent the intracolumn arcs with infinite cost, which guarantee that the desired surface intersects each column exactly once.The green inter-column arcs are constructed according to locally adaptive smoothness constraints learned from previous segmented surface.The yellow arcs are hard constraints with infinite cost.During optimization, the segmented surface will not cut across any hard constraint arcs.Therefore, the desired surface will not traverse the optic disc region, but is forced to go through the bottom of the optic disc region.Similarly, the hard constraint arcs are also constructed in y direction.For the multiple shared-hole surfaces segmentation, each graph is constructed for each surface in the same way as for single sharedhole surface graph search, and the inter-surface arcs are enforced to model the pairwise relations between surfaces.It is worth to mention that all the surface smoothness constraints and surface separation constraints are enforced outside the detected optic disc region.On the boundary of the optic disc region, only the hard constraints arcs are built to enforce the shared hole structure.Finally, the max-flow/min-cut algorithm is applied for optimization.
In this paper, we propose a novel shared-hole graph search method for ONH centered 3D SD-OCT image segmentation.Locally adaptive constraints are obtained from previously detected surfaces, thus avoiding the learning process, the detail of locally adaptive constraints are described in next section.

Image denoising
The dominating speckle noise in OCT scans may affect the effectiveness and efficiency of the following image analysis.However, most traditional denoising methods tend to blur away the boundaries in OCT images, and thus lead layer segmentation to trouble.In this work, anisotropic diffusion filter [42] is used to denoise the OCT volumes, which can reduce noise while preserving the edge information.

Surface 1 detection and B-scan alignment
The height of surface 1 changes drastically in ONH centered OCT image, causing trouble for conventional graph search method.In this paper, we introduce a novel two-stage method for surface 1 detection.The Otsu thresholding [43] combined with morphological operations is used for initialization, from which the search range and locally adaptive smoothness constraints are obtained.Then graph search is applied to refine the position of surface 1.
Otsu thresholding is an unsupervised method with automatic threshold selection for image segmentation.The method maximizes the intensity variance between the two classes, while minimizing the variance within each class at the same time.Morphological opening and closing procedures are applied to the Otsu segmented result to reduce noise and refine the segmented result.Finally, the upper boundary of the segmented result is regarded as the initialization of surface 1.The purpose of this step is to generate the shape and location constraints which are used to guide the graph search for refining surface 1.
The searching range are defined by the up-most constraint surface and down-most constraint surface, obtained as follows, ( ) ( ) where ( ) )  The eye movement during imaging often leads to discontinuity of retina structure in the 3D volume, causing difficulty for 3D analysis methods.As shown in Fig. 7 (a), the distortion can be perceived in y-z image, as the drastic jumping in adjacent columns corresponding to B-scans.To keep the natural curvature of retinal layers, we apply B-scan alignment [44] instead of image flattening [35,41,45].The average z position of the left most and right most 20% points in surface 1 in each B-scan is used to estimate the displacement of each B-scan.Each B-scan is thus shifted up or down to make sure that the average z positions of peripheral surface 1 become the same.

Polar-transformation of the volumetric image
In the original B-scans, the optic disc boundary is difficult to detect due to the presence of vessel shadows, especially for the slices where the optic disc boundaries get closer.Therefore, we polar-transform the Cartesian B-scans to radial B-scans, where the optic disc boundaries are kept apart, as the optic disc region are approximately round.The impact of vessel shadows thus becomes smaller and can be corrected in post-processing.To generate radial Bscans, the B-scan aligned volumetric image (Fig. 8(a)) is first interpolated in y direction into the same size as in x direction, making the x-y plane isotropic.Then, the center of x-y plane is regarded as the origin.The polar-transformation is applied inside the black circle (Fig. 8(b)), which ensures the generated radial B-scans are consistent in size.Finally, 180 radial B-scans, 1 degree apart, are generated with 2D linear interpolation.

Random forests for optic disc boundary detection
Random forests classifier [32] is an ensemble learning method for targets classification.The number of trees in the forest and the number of features randomly selected at each decision split are the key parameters for building a robust random forests classifier.Final classification is determined by taking the majority vote over the entire forests.In this paper, the number of trees and the number of features randomly selected at each decision split are set as 500 and 6.We extract 41 features (listed in Table 1) from polar-transformed image and classify each column in the polar-transformed image as inside or outside the optic disc region.Figure 9 shows a thumbnail of 41 extracted feature maps from a test data.The first feature is the intensity of projection image, which is produced by averaging the voxel intensities between surface 6 (IS/OS junction) and 9 (the outer boundary of the RPE).This project image contains the information of vessels and also indicates the location of the optic disc.Feature 2 is the average curvature of the surfaces 6, 8 and 9.The curvatures of intra-retinal surface indicate the structure change in the polar-transformed image.Feature 3 and 4 represent the difference before and after each voxel on the boundary of surface 9 in radial B-scan.Feature 3 is computed as the total intensity difference between two 30 × 30 patches to the left and right side of each voxel along surface 9 in the radial B-scan.As shown in Fig. 10(a), the purple line represents the detected boundary of surface 9 in radial B-scan, and the red dot is one voxel on surface 9.The feature 3 is the total voxels intensity value difference of yellow patch and green patch between right and left side of the red voxel.The value of this feature is much bigger in the boundary of optic disc region; as the intensity inside optic disc region (vitreous) is much lower than outside.Feature 4 is calculated as the total gradient difference 20 voxels before and after each voxel on surface 9.As shown in Fig. 10(b), the yellow and green segments are used to calculate the total gradient difference of the red voxel.The value of this feature is also much bigger in the boundary of optic disc region, as inside the optic disc region there is no layer structure and the gradient is small.Succinctly, feature 3 and 4 are calculated as follows,  G m B m are the intensity value and gradient value of the corresponding point in the radial B-scan, and M is the number of A-lines in the radial B-scan.Feature 5-6 are the location of the column, which deemed to be spatial prior information.Feature 7-9 are the height of surface 6, 8 and 9. Feature 10-41 are texture features extracted by Oriented Brief (ORB) algorithm [46], which is combines oriented scale-invariant feature transform (SIFT) [47] and rotation-aware brief (rBRIFE) [48].
Ideally, the boundary of ONH can be detected by solely using feature 4 [30].However, due to the strong noise of test data and the changeable ONH structure among glaucoma and normal eyes, feature 4 is not sufficient to detect the optic disc boundary accurately.Therefore, in this paper we introduce random forest classifier with the 41 features for accurately optic disc boundary detection.Some comparative experiments which indicate the effectiveness of the 41 features will be shown in the results section.
Refining process is needed to remove the influence of vessel shadows and noise.Morphological opening operator is first applied to eliminate small isolated points on the binary classification result.Then, a 2D graph search method [25] is used to detect the optic disc boundary based on the binary image.The smoothness constraint for 2D graph search is set to 1 to recover the smooth boundary.Then, the detection result is transformed back to original Cartesian B-scans. Figure 11 shows the optic disc boundary detection results both in the original and the radial B-scans.

Surface segmentation
In the proposed method, the multi-resolution strategy is used to improve the efficiency of surface segmentation.Figure 12 illustrates the down-sampling process.The workflows of surface segmentation are summarized as follows.9. Interpolate Surface 1-9 to original resolution data with smoothing.
Note that for surfaces 2-9, the search regions are limited by previously detected surfaces immediately above and below.For all surfaces, the parameters for computing adaptive smoothness constraints are set as 1

Data
In our study, the training data set for random forest included 3 ONH centered SD-OCT scans from normal eyes and 3 scans from glaucoma patients, and the test data set included 30 ONH centered SD-OCT scans from normal eyes and 35 scans from glaucoma patients.All OCT images were acquired using a Topcon 3D-OCT 2000 scanner (Topcon Corporation, Tokyo, Japan).Each volumetric image represented a 6mm × 6mm × 2.3mm region which contains 512 × 128 × 885 voxels ( X Y Z × × ).

Reference standard
For optic disc boundary detection, the reference standard came from one expert who annotated the optic disc boundary in each B-scan.The region inside the 2D projection of optic disc boundary was used as reference standard of optic disc region.Accuracy analysis was performed for all 65 test images.
For layer segmentation, from the test data sets, we randomly selected 10 OCT volume of normal eyes and 10 OCT volume of glaucomatous eyes.10 B-scans of each OCT volume which were centered at the central B-scan with interval of 5 B-scans were selected for retinal layer segmentation evaluation (200 scans in total).Two ophthalmologists independently traced the layer boundaries in each B-scan.The averages of these tracing were used as the reference standard.

Evaluation metrics and methods for comparison
To evaluate the optic disc region detection result, Dice similarity coefficient (DSC) was computed in the en-face (x-y) image to assess the accuracy of optic disc segmentation.The DSC measures the spatial overlap between two regions, A and B , defined as: We compared results of the proposed optic disc boundary segmentation method with Hu's method [25], Zang's method [30] and random forest trained with difference features.
To evaluate the surface segmentation results, the unsigned border positioning errors were calculated for each surface by measuring absolute Euclidean distances in the z-axis between segmentation results and the reference standard.For Surface 1 (ILM) segmentation, we compared the proposed method with the initialization method (Otsu thresholding combined with morphological operations), conventional multi-resolution graph search method with constant constraints [44], Zang's method [30], the Iowa Reference Algorithm by OCTExplorer software [53] (version 3.8.0),and GS-AAM method [31].Graph search method with constant constraints was applied on five resolution levels with down-sampled the OCT scan by a factor of 2 four times in z-direction, and the constant smoothness constraints for each levels on x-direction and y-direction were optimized empirically, set as 1, 1; 2, 2; 4, 4; 8, 8; 16, 16 from lowest resolution to original resolution.
For Surface 2-9 segmentation, we compared the proposed layer segmentation method with the conventional multi-resolution graph search method with constant constraints [44], Zang's method [30], the Iowa Reference Algorithm by OCTExplorer software [53] (version 3.8.0),and GS-AAM method [31].The multi-resolution segmentation strategy was set the same as proposed method.For surface 9, the smoothness constraints on x-direction and y-direction were set as 2, 2; 4, 4; 5, 5; 8, 8 from L5 Image to L1 Image.For surface 2-8, the smoothness constraints were set as 5, 5 and 8, 8 on L3 and L1 Image.The maximum and minimum surface separation constraints were set as 7, 4 for Surface 6 and 8; 25, 5 for Surface 2 and 3.These parameters were optimized empirically for the test data.
Note that in the comparison process, the surface 1 was compared in the entire range of Bscan, including the ONH region, and surface 2-9 were compared in the region with ONH excluded.For manual segmentation, the excluded ONH regions were based on the ground truth, while for automatic segmentation, the excluded ONH regions were based on the previous detected optic disc boundary.When there was a difference in the ONH boundaries between ground truth and detected results, the detected surfaces were linearly extrapolated to fill the missing part, or the redundant parts of detected surfaces were discarded.
Paired t-tests were performed both for optic disc detection and retinal surface segmentation, and a p-value less than 0.05 was considered statistically significant.

Optic disc region detection
To verify the effectiveness of the proposed 41 features, we ranked these features with permutation accuracy importance, which directly measured the impact of each feature on accuracy of the model [49].The general idea was to permute the values of each feature and measure how much the permutation decreases the accuracy of the model.Clearly, for unimportant features, the permutation should have little to no effect on model accuracy, while permuting important features should significantly decrease it.Figure 13 shows the permutation accuracy importance for each feature, the greater the permutation accuracy importance value, the more important of this feature is.With the rank of feature permutation accuracy importance, we compared proposed method with several different feature settings and several other methods: method 1, random forest trained with feature 1-9 (hand make features); method 2, random forest trained with feature 10-41 (ORB features); method 3, random forest trained with top 30 features of permutation accuracy importance; method 4, random forest trained with top 35 features of permutation accuracy importance; method 5, Zang's method [30] with feature 4; and method 6, Hu's method [25].Table 2 shows the DSC values of different methods.Based on the DSC values for the optic disc region detection, our classification based method outperformed Hu's method with statistically significant difference (p<0.001).Figure 14 shows examples of optic disc region detection results using proposed method and Hu's method.The first two rows show the detected disc region from glaucoma and normal data respectively, overlaid on the en-face projection images obtained by averaging the voxel intensities between surface 6 and 9.The third row shows the segmentation results in the central B-scan from the glaucoma data.The comparison demonstrates that the presence of the external oblique border tissue attaches to the end of BM surface causes erroneous optic disc boundary detection for Hu's method.However, our proposed method could handle this condition better.Compared to random forest trained with different features, the proposed method performed the best among different feature settings.As indicated by the results, ORB features are more important than the first 9 features, but the inclusion of these features improved the detection accuracy.Unfortunately, Zang's method [30] failed to detect optic disc boundary accurately on test data, probably due to the great noise in some test data and the changeable ONH structure among glaucoma and normal eyes.

Retinal surface segmentation
Table 3 shows the results of ILM segmentation with the unsigned border position error of the initialization by Otsu thresholding followed by morphological operations, the conventional multi-resolution graph search method [44], Zang's method [30], OCTExplorer software [53], GS-AAM method [31], and proposed method.Some segmentation results are shown in Fig. 15.The conventional multi-resolution graph search method [44] and GS-AAM method [31] failed to detect the steep slope inside the optic disc region due to constant smoothness constraint, resulting in large errors.Zang's method also failed to detect the boundary inside the optic disc region, because this method only considers 8 surrounding neighbors when determining the next optimal connection.The segmentation results of OCTExplorer software [53] were too smooth to accurately detect the ILM boundary inside the ONH region.The proposed initialization method resulted in a relatively good result, but was sometimes distracted by bright structures or artifacts inside the vitreous.The proposed Otsu segmentation guided graph search method learned the surface smoothness constraint from the initialization, refined it by global minimization of the gradient-based cost, and resulted in the best segmentation result both on steep and flat region.Table 4 shows the mean and standard deviation of unsigned border positioning errors for Surface 1-9, and the comparison results as p-values, where bold numbers indicate that the proposed method has statistically significantly better performance.The unsigned border positioning errors for 9 surfaces segmentation indicated that our proposed method performed better than other 4 compared methods.
Surface segmentation examples obtained by conventional multi-resolution graph search method [44], Zang's method [30], OCTExplorer software [52], GS-AAM method [31] and the proposed method are shown in Fig. 15.The conventional multi-resolution graph search segmentation results ran arbitrary inside the ONH region, which caused big errors both on the adjacent columns of the same B-scan and on the adjacent B-scans, due to the usage of constant smoothness constraints in both x and y-direction.On the other hand, the detected surface 1 sometimes fell to surface 6, due to the usage of large constant smoothness constants.Zang's method only considered 8 surrounding neighbors when determining the next optimal connection, which also failed to detect the steep boundary.And this 2D method detected retinal layers slice by slice, and no multiple layer segmentation strategy was applied, these shortcomings made the segmentation layers interfered by nearby undetected layers and resulted in error segmentation.OCTExplorer software [53] showed smooth surface curves, but the biggest problem was OCTExplorer software [53] could not accurately segment the surfaces near the deep structure region.This is because the deep structure inside the ONH region misleads the segmented surfaces away from their actual locations.And the GS-AAM method [31] is unable to learn the changeable ONH structure among glaucoma and normal eyes, which produced wrong information for graph search method, and failed to detected accurate surfaces.However, with the proposed shared-hole graph construction and locally adaptive constraints, the proposed method could achieve better results.Besides optic disc region, the optic cup region also plays an important role in the glaucoma diagnosis.In our subsequent studies, we will add optic cup analysis, so that the cup-to-disc ratio can be calculated for quantitative analysis of glaucoma.

Conclusion
In summary, we have proposed a novel framework for optic disc region detection and 9 retinal surfaces segmentation on ONH centered SD-OCT scans.Our layer segmentation method achieved accurate segmentation by improving the graph search algorithm with shared-hole and locally adaptive constraints.Our proposed random forest classifier trained with novel surface-derived features achieved higher optic disc region detection accuracy compared to the state-of-the-art methods.A mean Dice coefficient of 0.925 ± 0.03 was achieved for optic disc region detection, and an overall mean unsigned border positioning error of 7.27 ± 5.40 µm was achieved for 9 surfaces segmentation.The experimental results showed the feasibility of applying optic disc boundary detection and retinal layers segmentation in clinical diagnosis.

Fig. 4 .
Fig. 4. Graph construction for graph search with surface smoothness constraints (a) and surface separation constraints (b).
(b) illustrates the surface separation constraints for a pair of surfaces.The maximum surface separation constraint is set as

Fig. 5 .
Fig. 5. Graph construction for single surface shared-hole graph search in the x-z plane.The translucent yellow region indicates the shared-hole region for graph search.The yellow arcs are the hard constraint arcs.

2 x d = , 1 yd 4 x
adaptive smoothness constraints for each pair of neighboring columns in the x and y-direction respectively.x d and y d are the default smoothness constraints for graph search (set to = empirically), which prevent the smoothness constraints becoming too small.Figure6shows an example of the results obtained by the proposed method, Otsu thresholding, and conventional graph search method with small constant smoothness constraints ( .The result of small smoothness constraints fails to follow the boundary in the deep disc cup region (Fig.6(a)), and the result of big smoothness constraints is easily attracted by noise (Fig.6(c)).The initialization based on Otsu thresholding performs better, but deviates a bit when there are high intensity structures above surface 1(Fig.6(b)).The proposed method detects surface 1 accurately in all locations.

Fig. 6 .Fig. 7 .
Fig. 6.Surface 1 segmentation results.The yellow and pink curves are the results of conventional graph search with small and big constant constraints, the red curve is the initialization result by Otsu thresholding with morphological operation, and the green curve shows the result of Otsu segmentation guided graph search with locally adaptive constraints.
Figure 7(b) shows the result after B-scan alignment.

Fig. 8 .
Fig. 8.The polar-transformation process.B-scan aligned volumetric image (a) is first interpolated in y direction into the same size as in x direction, making the x-y plane isotropic.Then, the center of x-y plane is regarded as the origin.Finally, inside the black circle area (b), 180 radial B-scans, 1 degree apart, is generated with 2D linear interpolation.(c) Three polar transformed radial B-scans on different angles.

Fig. 10 .
Fig. 10.Sketch map for feature 3 and 4 in a radial B-scan.(a) Feature 3 calculated the total intensity difference between two 30 × 30 patches (yellow patch and green patch) to the left and right side of each voxel (red dot) along surface 9 (purple curve) in the radial B-scan.(b) Feature 4 is calculated as the total gradient difference 20 voxels (yellow curve and green curve) before and after each voxel (red dot) on surface 9 (purple curve).
the position of detected boundary of surface 9 in A-line i and m , ( , ) I m n and ( , ( ))

Fig. 11 .
Fig. 11.Optic disc boundary detection results.Left: results showed on original B-scans, Right: results showed on radial B-scans.

1 . 4 . 6 .
Detect Surface 1 in L1 Image by Otsu thresholding guided graph search with locally adaptive constraints.2. Detect Surface 9 in L5 Image by shared-hole single surface graph search, with locally adaptive constraints learned from down-sampled Surface 1.Then, refine Surface 9 in L4 Image by shared-hole single surface graph search, with locally adaptive constraints learned from interpolated initial detection of Surface 9.This refinement is repeated in L2 and L1 Image consecutively. 3. Detect Surface 6 and 8 in L3 Image by shared-hole double surface graph search, with locally adaptive constraints learned from detected Surface 9 in L3 Image.The maximum surface separation constraint is set as max , Detect Surface 4 in L3 Image by shared-hole single surface graph search, with locally adaptive constraints learned from down-sampled Surface 1 in L3 Image. 5. Detect Surface 2 and 3 in L3 Image by shared-hole double surface graph search, with locally adaptive constraints learned from detected Surface 4 in L3 Image, and The maximum surface separation constraint is set as max , Detect Surface 5 in L3 Image by shared-hole single surface graph search, with locally adaptive constraints learned from detected Surface 4 in L3 Image.7. Detect Surface 7 in L3 Image by shared-hole single surface graph search, with locally adaptive constraints learned from detected Surface 9 in L3 Image.8. Refine Surface 2-9 in L1 Image by shared-hole single surface graph search, with locally adaptive constraints learned from interpolated Surface 2-9 from L3 Image respectively.

Fig. 14 .
Fig. 14.An example of optic disc detection result.The first two rows show the detected disc region from glaucoma and normal data respectively, (a) (e) the reference standard, (b) (f) the result of the proposed method, (c) (g) the result of Hu's method, (d) (h) the boundary overlay of reference standard, proposed method and Hu's method.The third row shows a central Bscan from the glaucoma data as used in first row.(i) the reference standard, (j) the result of the proposed method, (k) the result of Hu's method [25].

Table 4 . Summary of Mean Unsigned Border Position Error for All Data
Mean ±SD in μm, 2.6μm = 1 pixel.