Automated segmentation of optic disc in SD-OCT images and cup-to-disc ratios quantification by patch searching-based neural canal opening detection

Glaucoma is one of the most common causes of blindness worldwide. Early detection of glaucoma is traditionally based on assessment of the cup-to-disc (C/D) ratio, an important indicator of structural changes to the optic nerve head. Here, we present an automated optic disc segmentation algorithm in 3-D spectral domain optical coherence tomography (SD-OCT) volumes to quantify this ratio. The proposed algorithm utilizes a two-stage strategy. First, it detects the neural canal opening (NCO) by finding the points with maximum curvature on the retinal pigment epithelium (RPE) boundary with a spatial correlation smoothness constraint on consecutive B-scans, and it approximately locates the coarse disc margin in the projection image using convex hull fitting. Then, a patch searching procedure using a probabilistic support vector machine (SVM) classifier finds the most likely patch with the NCO in its center in order to refine the segmentation result. Thus, a reference plane can be determined to calculate the C/D radio. Experimental results on 42 SDOCT volumes from 17 glaucoma patients demonstrate that the proposed algorithm can achieve high segmentation accuracy and a low C/D ratio evaluation error. The unsigned border error for optic disc segmentation and the evaluation error for C/D ratio comparing with manual segmentation are 2.216 ± 1.406 pixels (0.067 ± 0.042 mm) and 0.045 ± 0.033, respectively. ©2015 Optical Society of America OCIS codes: (100.0100) Image processing; (110.4500) Optical coherence tomography; (170.4470) Ophthalmology. References and links 1. N. G. Strouthidis, H. Yang, B. Fortune, J. C. Downs, and C. F. Burgoyne, “Detection of optic nerve head neural canal opening within histomorphometric and spectral domain optical coherence tomography data sets,” Invest. Ophthalmol. Vis. Sci. 50(1), 214–223 (2009). 2. J. S. Schuman, G. Wollstein, T. Farra, E. Hertzmark, A. Aydin, J. G. Fujimoto, and L. A. Paunescu, “Comparison of optic nerve head measurements obtained by optical coherence tomography and confocal scanning laser ophthalmoscopy,” Am. J. Ophthalmol. 135(4), 504–512 (2003). 3. Y. H. Kwon, M. Adix, M. B. Zimmerman, S. Piette, E. C. Greenlee, W. L. Alward, and M. D. Abràmoff, “Variance owing to observer, repeat imaging, and fundus camera type on cup-to-disc ratio estimates by stereo planimetry,” J. Glaucoma 18(4), 305–310 (2009). 4. M. D. Abràmoff, K. Lee, M. Niemeijer, W. L. Alward, E. C. Greenlee, M. K. Garvin, M. Sonka, and Y. H. Kwon, “Automated segmentation of the cup and rim from spectral domain OCT of the optic nerve head,” Invest. Ophthalmol. Vis. Sci. 50(12), 5778–5784 (2009). #251419 Received 5 Oct 2015; revised 16 Nov 2015; accepted 17 Nov 2015; published 20 Nov 2015 © 2015 OSA 30 Nov 2015 | Vol. 23, No. 24 | DOI:10.1364/OE.23.031216 | OPTICS EXPRESS 31216 5. J. Xu, H. Ishikawa, G. Wollstein, R. A. Bilonick, L. Kagemann, J. E. Craig, D. A. Mackey, A. W. Hewitt, and J. S. Schuman, “Automated volumetric evaluation of stereoscopic disc photography,” Opt. Express 18(11), 11347– 11359 (2010). 6. A. Aquino, M. E. Gegúndez-Arias, and D. Marín, “Detecting the optic disc boundary in digital fundus images using morphological, edge detection, and feature extraction techniques,” IEEE Trans. Med. Imaging 29(11), 1860–1869 (2010). 7. H. Yu, E. S. Barriga, C. Agurto, S. Echegaray, M. S. Pattichis, W. Bauman, and P. Soliz, “Fast localization and segmentation of optic disk in retinal images using directional matched filtering and level sets,” IEEE Trans. Inf. Technol. Biomed. 16(4), 644–657 (2012). 8. J. Cheng, J. Liu, Y. Xu, F. Yin, D. W. K. Wong, N. M. Tan, D. Tao, C. Y. Cheng, T. Aung, and T. Y. Wong, “Superpixel classification based optic disc and optic cup segmentation for glaucoma screening,” IEEE Trans. Med. Imaging 32(6), 1019–1032 (2013). 9. C. Muramatsu, T. Nakagawa, A. Sawada, Y. Hatanaka, T. Hara, T. Yamamoto, and H. Fujita, “Automated segmentation of optic disc region on retinal fundus photographs: Comparison of contour modeling and pixel classification methods,” Comput. Methods Programs Biomed. 101(1), 23–32 (2011). 10. 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). 11. J. Xu, H. Ishikawa, G. Wollstein, and J. S. Schuman, “3D optical coherence tomography super pixel with machine classifier analysis for glaucoma detection,” in Proceedings of IEEE conference on Engineering in Medicine and Biology Society (IEEE, 2011), pp. 3395–3398. 12. 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 Medical Imaging 72620U (2009). 13. K. Lee, M. Niemeijer, M. K. Garvin, Y. H. Kwon, M. Sonka, and M. D. Abràmoff, “3-D segmentation of the rim and cup in spectral-domain optical coherence tomography volumes of the optic nerve head,” Proc. SPIE 7262, 72622D (2009). 14. K. Lee, M. Niemeijer, M. K. Garvin, Y. H. Kwon, M. Sonka, and M. D. Abràmoff, “Segmentation of the optic disc in 3-D OCT scans of the optic nerve head,” IEEE Trans. Med. Imaging 29(1), 159–168 (2010). 15. M. S. Miri, K. Lee, M. Niemeijer, M. D. Abràmoff, Y. H. Kwon, and M. K. Garvin, “Multimodal segmentation of optic disc and cup from stereo fundus and SD-OCT images,” Proc. SPIE Medical Imaging 86690O (2013). 16. N. G. Strouthidis, H. Yang, J. F. Reynaud, J. L. Grimm, S. K. Gardiner, B. Fortune, and C. F. Burgoyne, “Comparison of clinical and spectral domain optical coherence tomography optic disc margin anatomy,” Invest. Ophthalmol. Vis. Sci. 50(10), 4709–4718 (2009). 17. Z. Hu, M. D. Abràmoff, Y. H. Kwon, K. Lee, and M. K. Garvin, “Automated segmentation of neural canal opening and optic cup in 3D spectral optical coherence tomography volumes of the optic nerve head,” Invest. Ophthalmol. Vis. Sci. 51(11), 5708–5717 (2010). 18. F. A. Medeiros, L. M. Zangwill, C. Bowd, R. M. Vessani, R. Susanna, Jr., and R. N. Weinreb, “Evaluation of retinal nerve fiber layer, optic nerve head, and macular thickness measurements for glaucoma detection using optical coherence tomography,” Am. J. Ophthalmol. 139(1), 44–55 (2005). 19. A. Manassakorn, K. Nouri-Mahdavi, and J. Caprioli, “Comparison of retinal nerve fiber layer thickness and optic disk algorithms with optical coherence tomography to detect glaucoma,” Am. J. Ophthalmol. 141(1), 105–115 (2006). 20. T. C. Chen, “Spectral domain optical coherence tomography in glaucoma: qualitative and quantitative analysis of the optic nerve head and retinal nerve fiber layer (an AOS thesis),” Trans. Am. Ophthalmol. Soc. 107, 254–281 (2009). 21. 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). 22. P. Hrynchak, N. Hutchings, D. Jones, and T. Simpson, “A comparison of cup-to-disc ratio measurement in normal subjects using optical coherence tomography image analysis of the optic nerve head and stereo fundus biomicroscopy,” Ophthalmic Physiol. Opt. 24(6), 543–550 (2004). 23. Q. Chen, T. Leng, L. Zheng, L. Kutzscher, J. Ma, L. de Sisternes, and D. L. Rubin, “Automated drusen segmentation and quantification in SD-OCT images,” Med. Image Anal. 17(8), 1058–1072 (2013). 24. M. Gargesha, M. W. Jenkins, A. M. Rollins, and D. L. Wilson, “Denoising and 4D visualization of OCT images,” Opt. Express 16(16), 12313–12333 (2008). 25. C. Tomasi and R. Manduchi, “Bilateral filtering for gray and color images,” in Proceedings of IEEE Conference on Computer Vision (IEEE, 1998), pp. 839–846. 26. Q. Chen, W. Fan, S. Niu, J. Shi, H. Shen, and S. Yuan, “Automated choroid segmentation based on gradual intensity distance in HD-OCT images,” Opt. Express 23(7), 8974–8994 (2015). 27. 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). 28. 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). 29. X. Chen, J. K. Udupa, U. Bagci, Y. Zhuge, and J. Yao, “Medical image segmentation by combining graph cuts #251419 Received 5 Oct 2015; revised 16 Nov 2015; accepted 17 Nov 2015; published 20 Nov 2015 © 2015 OSA 30 Nov 2015 | Vol. 23, No. 24 | DOI:10.1364/OE.23.031216 | OPTICS EXPRESS 31217 and oriented active appearance models,” IEEE Trans. Image Process. 21(4), 2035–2046 (2012). 30. C. Cortes and V. Vapnik, “Support-vector networks,” Mach. Learn. 20(3), 273–297 (1995). 31. T. Ojala, M. Pietikäinen, and T. Mäenpää, “Multiresolution gray-scale and rotation invariant texture classification with local binary patterns,” IEEE Trans. Pattern Anal. 24(7), 971–987 (2002). 32. T. Ojala, M. Pietikäinen, and T. Mäenpää, “A generalized Local Binary Pattern operator for multiresolution gray scale and rotation invariant texture classification,” in Proceedings of International Conference on Advances in Pattern Recognition (IEEE, 2001), pp. 399–406. 33. N. Dalal and B. Triggs, “Histograms of oriented gradients for human detection,” in proceedings of IEEE Conference on Computer Vision and Pattern Recognition (IEEE, 2005), pp. 886–893. 34. J. Zhang, K. Huang, Y. Yu, and T. Tan, “Boosted local structured hog-lbp for object localization,” in proceedings of IEEE Conference on Computer Vision and Pattern Recognition (IEEE, 2011), pp. 1393–1400. 35. C. C. Chang and C. J. Lin, “LIBSVM: A library for support vector machines,” ACM. Trans. Intel. Syst. Technol. 2(3), 1–27 (2011). 36. G. J. Jaffe and J. Caprioli, “Optical coherence tomography to detect and manage retinal disease and gl


Introduction
Glaucoma is a chronic neurodegenerative disease of the optic nerve, which demonstrates classic structural characteristics including cupping of the optic nerve head, focal and diffuse retinal rim loss, and nerve fiber layer defects [1].There is a need to identify early internal structural changes in glaucoma.The cup-to-disc (C/D) ratio is an important indicator for evaluating the glaucomatous changes of the optic nerve head (ONH) [2].To quantify the ratio, specialists commonly perform planimetry in 2-D color fundus photographs.However, manual planimetry has been proven time-consuming and is wrought with interobserver variability [3,4].With the introduction of spectral-domain optical coherence tomography (SD-OCT), which is a high-resolution cross-section and noncontact imaging technology, it is possible to image the ONH in a 3-D manner, enabling a more detailed quantification of its structure and possible changes [5].
In order to quantify the C/D ratio, many automatic methods have segmented the optic disc and cup in color fundus photographs.For example, Aquino et al. proposed a template-based method to detect optic disc boundary, using circular Hough transform for boundary approximation [6].Yu et al. also presented a hybrid level-set approach for optic disc segmentation based on the deformable model, which combined region and local gradient information [7].In [8], a superpixel classification based algorithm was applied to determine the optic disc by contrast enhanced histogram and center surround statistics.Furthermore, comparisons of the active contour models for glaucoma screening were performed in [9].
Recently, there have been several studies for glaucoma detection in SD-OCT images [10][11][12][13][14][15].Antony et al. presented an automated intraretinal layer segmentation algorithm to calculate the thickness of the retinal nerve fiber layer in normal and glaucoma scans [10].To detect glaucoma structural damage at an early stage, Xu et al. generated a 2-D feature map from a SD-OCT volume by grouping super pixels and utilized a boosting algorithm to classify glaucoma cases [11].Another category of approaches is based on optic disc segmentation.Work by Hu et al. [12] transformed the SD-OCT slices to planar projection images, and used graph search algorithm to detect the two boundaries of optic disc and cup simultaneously.In [13,14], Lee et al. proposed a multi-scale 3-D graph search algorithm to segment retinal surfaces for OCT scan flattening, and then classified each voxel column (A-scan) using k-NN classifier according to the features obtained from the projection image and the retinal surfaces.Based on this work, Miri et al. proposed a multi-modal pixel classification method to segment the optic disc, combining stereo fundus and SD-OCT volumes [15].However, these previous methods required the assistance of color fundus photographs.Additionally, A-scanbased classification may suffer from morphological diversity of the ONH and lose important 3-D global information.
In this paper, we propose an automated optic disc segmentation and C/D ratio quantification method based on neural canal opening (NCO) detection in the optic nerve head.In [1,16], Strouthidis et al. pointed out that NCO is an objective anatomic landmark consistent with SD-OCT optic disc margin anatomy.Furthermore, Hu et al. found out that the NCO can serve as longitudinally stable reference plane and it is not likely to change with glaucomatous progression [17].As shown in Fig. 1(a), NCO is defined as the termination of the retinal pigment epithelium (RPE) layer.A parallel line of 150 μm above the connecting line of the two NCO points indicates the standard reference plane [18,19].The circular ring between two intersections of the ONH surface and the reference plane are defined as the cup border.Although there has been some debate about the determination of the reference plane [20][21][22], we adopted this classic definition.In Fig. 2(b), the red and the green dots are the points of disc and cup margin in projection fundus image, corresponding to the columns of the NCO and the cup borders, respectively.Therefore, the disc margin can be constituted by the NCO in SD-OCT images.
Unlike the methods that directly segment the optic disc in projection image or using classifier to determine which A-scan belongs to the cup or rim, our approach attempts to extract the NCO from SD-OCT scans for the delineation of the disc margin.The proposed approach utilizes a two-stage strategy.The first step is to locate the coarse disc margin by the segmentation of the retinal pigment epithelium (RPE) layer and the smooth constraint of consecutive B-scans.In the second step, we develop a support vector machine (SVM)-based patch search method to find the most likely patch centered at the NCO and refine the segmentation result.Using the NCO and reference plane, the cup border can be evaluated.Finally, the C/D ratio is calculated by the cup diameter dividing the disc diameter.To the best of our knowledge, our approach is the first to automatically segment the optic disc by a NCO detection-based patching search method.The two-stage strategy combines the global and local information for optic disc segmentation.When determining the coarse disc margin location, we utilize the structural characteristics of ONH for initial NCO detection, while in the patch searching procedure, the visual content similarity near the NCO is applied for the final segmentation.

Overview of the method
Figure 2 shows the flowchart of the proposed algorithm, which comprises two main stages: the coarse disc margin location and the SVM-based patch search.In the first stage, each B-scan in the volume is denoised and rescaled during preprocessing.Then, a 3-D graph search algorithm is applied to automatically segment the RPE layer [14].Based on the segmentation result, we determined the initial NCO position by the maximum curvature of the detected RPE boundary and the smooth spatial constraint of the consecutive B-scans.In the second stage, we select image patches from SD-OCT volumes, and utilize a probabilistic SVM classifier for training after feature extraction of the patches.Then, the searching procedure is generated at the region restricted by the initial NCO location.The patch of maximum probability, centered at NCO, is regarded as the final NCO position.After the two steps, the cup border can be calculated by the location of the NCO and the ILM boundary, and the C/D ratio can be quantified.

Coarse disc margin location
The purpose of this step is to limit NCO region for the subsequent patch search procedure.This step involves image preprocessing, RPE segmentation and initial NCO detection.

Preprocessing
Since coherent optical beams have low temporal coherence and high spatial coherence, SD-OCT images contain speckle noise [23].Speckle size may vary in the axial and lateral dimensions, which is mainly determined by source bandwidth and numerical aperture.In addition, shot noise is also present in SD-OCT images, which can be adequately described by the additive white Gaussian noise (AWGN) process [24].To facilitate the segmentation of RPE and ILM layer, we pre-processed the images by bilateral filtering denoising [25,26], which has proven adequate for edge preservation.
For added effectiveness in the NCO detection, we also resampled the SD-OCT images in x -axis direction.Specifically, each B-scan in the volume was rescaled from 200 × 1024 to 600 × 400 in the axial and x -axis direction, respectively.An example of a B-scan after denoising and rescaling is shown in Fig. 3(b).

RPE layer segmentation
The 3-D graph search-based segmentation algorithm has proven to be effective for detecting multiple intraretinal surfaces [27][28][29].With some appropriate feasibility constraints, the surface segmentation problem can be directly converted to a minimum closed set finding problem in a geometric graph and solved in a low-order polynomial time.We applied the multi-scale 3-D graph search approach introduced in [14] to segment the outer boundary of the RPE for subsequent NCO detection.As shown in Fig. 3(c), the green and yellow curves indicate the segmented ILM and RPE boundaries, respectively.

Initial NCO detection
From the RPE curve fitting in Fig. 3(c), we observed that the visible parts of the RPE follow an approximately linear curvature, and change dramatically at the termination of the RPE.Therefore, we selected points with maximum curvature as the NCO.However, there are several morphological characteristics that should be noticed: (1) The NCO points are restricted in a certain region by the imaging position and prior anatomical knowledge.The NCO always localizes near the center of the image, as the images are taken aimed to be centered at the optic disc, and it is also noticeable that the axial depth of the NCO locations do not vary substantially from the average height of the visible parts of RPE.(2) The position of the NCO locations should not appreciably change in consecutive B-scans.In this paper we refer to this as a spatial correlation smoothness constraint.Based on these constraints, we firstly calculated the positions of NCO candidates as: where 0 p denotes the NCO candidate in the RPE boundary and ( , ) C x y represents the curvature of the given point.We used the location constraint to ensure 0 p in an appropriate region.w denotes the width of the image so that the lateral position x is limited at the center part of the image.We also restricted the height of the NCO between 1 h and 2 h , which are two adaptive parameters.Taking the left NCO point as an example, we fitted a straight line to the left-side flat part of the RPE boundary, indicated as ( )  g x (shown in red in Fig. 3 , where 1 x and 2 x are assumed as the start and the end x -coordinate of the fitting line, fixed to 0 and 300 in this paper.The equivalent operation was also conducted for the right NCO point and the right-side part of the RPE boundary.An example of the initial NCO detection is shown in Fig. 3(d).
In the projection fundus image, we firstly defined a scope of disc margin in x -axis direction as [ 10,10] l r , where l x and r x are the average x -coordinate of the left and right NCO candidates in B-scans 95#~105#, respectively, considering the distance between the left and right NCO is largest in central B-scans, and filtered the individual outliers beyond this scope.Then, we defined a disc margin smoothing procedure by the spatial correlation smoothness constraint as: where i p indicates the point of the disc margin in projection image corresponding to the column of 0 p , and i y is the serial number of the B-scan.The ( ) N p denotes the disc margin points in former 2 and latter 2 consecutive B-scans of the current B-scan defined as a neighborhood and 4 n = represents the neighborhood size, while the constant k is set to 5. The x -coordinate locations of the initially selected points are smoothed by considering the selected locations in neighboring B-scans, in order to restrict the NCO locations in the consecutive B-scans to not change more than a threshold amount.
Figure 3(e) shows the NCO detection result in the projection fundus image.The green curve indicates the coarse disc margin by mapping NCO points from each B-scan to the projection image.For preserving the shape of the optic head, we refined the disc margin using a convex hull fitting operation, as shown in Fig. 3(f).A manual expert-defined reference standard is shown in Fig. 3(g) for comparison.
Although the result indicates that the coarse disc margin is close to the reference standard, its determination still depends on the precision of the RPE segmentation.The points with maximum curvature may not be the optimal location for the NCO when the segmentation result is poor.Therefore, we developed a SVM-based patch searching method to refine the NCO detection in a following step.

SVM-based patch searching
Because the initial NCO detection is influenced by the RPE segmentation, we developed a patch searching method to refine the optic disc segmentation.The purpose of this method is to find the most likely patch whose center represents the NCO.We utilized a SVM classifier [30] to select the patches with maximum probability of belonging to classes defined as those centered at the left or right NCO (Figs. 4(b) and 4(c)), obtained from slide windows near the initial NCO location (Fig. 4(a)).4(e)).The size of each patch was set to be 81 × 81 to ensure that it is discriminative and robust for classification.All the patches were selected near the true NCO or NCO candidates calculated in section 2.2.3.

Feature extraction
We extracted two kinds of texture features for patch description: local binary pattern (LBP), histogram of gradient (HOG).The LBP and HOG features were then combined to form a complete feature set.
The LBP is a simple and efficient textural operator which labels the pixels of an image by thresholding the neighborhood of each pixel with the value of the center pixel and generates the result as a binary number [31,32].It can be viewed as a unifying approach to the traditionally divergent statistical and structural models of texture analysis.In this paper, we applied a typical circle LBP operator introduced in [32] for computational efficiency, where the number of neighborhood and the radius are 8 and 2, respectively.Therefore, we obtained a 59 dimensional feature vector of LBP for each image sample.
The HOG descriptor counts the occurrences of gradient orientation in localized portions of an image [33,34].It considers that the local object appearance and shape within an image can be described by the distribution of intensity gradients or edge directions.Because of its discriminative power and computational simplicity, HOG has become a popular approach in various applications.An 81 dimensional HOG feature vector was calculated for each image sample here.
We used a serial fusion strategy to combine LBP and HOG features.By a serial linear combination, the two types of the feature vectors are fused into a discriminant vector for classification.Hence the dimension of the fusion feature is 140.

Patch searching
After SVM training, we determined a searching range and calculated the probability of each patch within the range belonging to class 1 or 2. We assumed that the true NCO is near the NCO candidate selected in the initial detection step, so we limited the search range in the xdirection to 0 0 [ 30,30] x x − + pixels, where 0 x denotes the x -coordinate of the NCO candidate.To reduce the influence of RPE segmentation error, we allowed a height offset of [-5, 10] pixels from the detected RPE boundary.We also defined an interval between two slide windows of 5 pixels.For achieving more precise segmentation results, we finely imposed a set of constrains and refinements in the optimal patch selection.Firstly, we assumed that the calculated NCO should be close enough to the initial NCO candidates.If the distance between these two points is more than 35 pixels, we select the top 5 patches with maximum probability of belonging to class 1 or 2, and refined NCO location as the center of the patch closest to the NCO candidate.Secondly, according to the spatial correlation smoothness constraint of the consecutive Bscans, we adjusted the position of NCO by Eq. ( 2).

Experiments and results
We acquired 42 SD-OCT volumes (20 left eye and 22 right eye) centered at optic nerve heads from 17 glaucoma patients to test our automated segmentation algorithm.All of the volumes were obtained using a Cirrus SD-OCT device (Carl Zeiss Meditec, Inc., Dublin, CA).Each volume consisted of 200 contiguous 200 × 1024 voxel B-scan images (each B-scan comprising 200 A-scans containing 1024 voxels).The voxel size was 30 × 30 × 2 μm and its depth was 8 bits in grayscale.
The proposed algorithm was implemented in Matlab 2010b and run on a 3.30GHz Intel Xeon PC with 16 GB memory.We selected 1600 independent samples from 20 SD-OCT volumes (400 samples for each class) as the training set.The libsvm software package with RBF kernel was applied to train the probabilistic SVM classifier [35].The penalty factor C and the kernel parameter γ were determined by cross validation for different kinds of features (LBP: C = 5.0, γ = 0.49; HOG: C = 5.0, γ = 0.07; Fusion: C = 5.0, γ = 0.30).

Evaluation of optic disc segmentation
To evaluate the performance of the proposed algorithm, we compared our segmentation results with a previously proposed A-scan classification-based segmentation method and a manual segmentation on the 42 SD-OCT test volumes.Segmentation methods based on Ascan classification have been previously proposed in [13][14][15] and aim to label each A-scan (corresponding to one pixel in the projection image) in B-scan images as cup, rim or background using k-NN classifier.We extracted a 15 dimensional feature vector from each Ascan as described in [14] to train the k-NN classifier.The manual segmentations were generated by two experienced experts who manually marked the NCO for each B-scan and calculated the cup border by the reference plane to segment optic disc and cup in projection images.The reference standard was obtained based on the two readers' consensus results in projection images.
We utilized unsigned border error (UBE) and Dice similarity coefficient (DSC) to estimate the accuracy of the tested segmentation algorithms [13,14].The UBE indicates the average closest distances between all boundary points from segmentation regions of an algorithm and reference standard, while the DSC denotes the spatial overlap between those two regions.Considering S and R as the regions outlined by a segmentation algorithms and a reference standard, respectively, the UBE and DSC were calculated as: ( ) where the min ( , ) Dist r s denotes the minimum Euclidean distance between one pixel in region r and all the pixels in region s , and ( ) Num a is the number of pixels in region a .Table 1 and Table 2 report the UBE and DSC of the disc and cup segmentations using the algorithm presented here and different feature sets (LBP, HOG, and fusion of both of them).It is observed that the performance using the fusion of LBP and HOG features is slightly better than LBP and HOG, which indicates that the fusion features are more discriminative for our segmentation algorithm.Therefore, we selected the fusion features for subsequent evaluation.As a qualitative evaluation, Fig. 5 displays the NCO detection by our algorithms in B-scan images and the comparisons of ONH segmentation by different algorithms in projection images.It was apparent that our patch searching-based method achieved a more accurate disc and cup segmentation than the coarse disc margin location method (introduced in section 2.2.3) and A-scan classification-based segmentation.The optic disc margin segmented by coarse disc margin location seems consistently larger than by the manual segmentation since it applies convex hull fitting to preserve the largest margin for the patch searching procedure.The disc and cup margins segmented by the patch searching-based method are smoother than by the A-scan classification-based segmentation because it considers the spatial correlation smoothness constraint, rather than by directly classifying each A-scan independently.Table 3 and Table 4 show the UBE and DSC of the different evaluated algorithms comparing with the reference standard from the 42 SD-OCT test volumes.

Run time comparison
To evaluate the efficiency of optic disc segmentation, we recorded the average computational time of our algorithm with different features and compared this with that from the A-scan classification-based segmentation algorithm, shown in Table 6.We determined that the run time of patch searching is mainly determined by feature extraction.Although fusion feature is more discriminating for our algorithm, it required more time to perform segmentation

Conclusion and discussion
In this paper, a patch search-based optic disc segmentation algorithm has been proposed for quantifying the C/D ratio in SD-OCT volumes.Compared with traditional methods, our algorithm has several advantages: (1) A two-stage strategy combines the global and local information for optic disc segmentation.In coarse disc margin location, we consider the structural characteristics of ONH for initial NCO detection.In the patch searching procedure, we seek to find the most likely patch centered at the NCO near NCO candidates.(2) Unlike A-scan classification-based segmentation methods that directly classify every column in Bscan images as disc (rim or background), our algorithm only searches in a restricted region, increasing efficiency.Furthermore, the SVM classifier spends less computational time than a k-NN classifier because it can use pre-training model data, further boosting efficiency.(3) As shown in Fig. 5, the segmentation results of the proposed algorithm have smoother margins than those of other methods.This is due to the fact that our method imposes spatial correlation smoothness constraints to finely tune the final segmented result.Although the experimental results demonstrate that the proposed algorithm can achieve high segmentation accuracy and be an efficient clinical tool for quantifying the C/D ratio, there are several limitations: (1) The computation cost principally depends on feature extraction.In Table 6, it is observed that the time cost of our algorithm with the fusion feature does not yet meet the needs of real-time segmentation.Therefore, efficient feature extraction is required.(2) The results produced by our method could present inaccuracies in severely tilted images.Figure 7 shows an example of inaccurate NCO detection.The patch has a maximum probability of class 2 in the searching procedure, but the detected NCO is far away from the true NCO since the similarity of visual content may not fully reflect anatomical characteristics.
In future research, we plan to improve our algorithm in two ways.Feature selection will be refined to increase the computational efficiency and NCO detection will be improved by incorporating anatomical structures to improve the segmentation accuracy.

Fig. 1 .
Fig. 1.(a) Central B-scan (number 100 of 200) of an example SD-OCT volume.The two NCO points at the end of the RPE (indicated by the yellow outline as automatically segmented by 3-D graph search algorithm) are denoted by red dots.The green line indicates the reference plane and its intersections with inner limiting membrane (ILM, indicated by the blue outline as automatically segmented by 3-D graph search algorithm) are cup boundary points denoted by green circle markers.(b) The projection fundus image of the same volume.The red and green dots are the points of optic disc and cup margin, respectively, which correspond to the columns of NCO and cup borders in (a).

Fig. 4 .
Fig. 4. SVM-based patch searching.(a) Patch searching using slide windows denoted by red dotted boxes.(b) -(e) four classes of the image samples.

2. 3 . 1
Sample selection For SVM classifier training, we selected 1600 image samples from 20 SD-OCT volumes as training set, which were divided into four classes: (1) patches centered at the left NCO (Fig. 4(b)); (2) patches centered at the right NCO (Fig. 4(c)); (3) patches including the RPE layer (Fig. 4(d)); (4) patches of background among two NCO points (Fig. vectors obtained from the patches, the probability that a patch j I belong to class k is determined by the prediction of SVM as jk p .The maximum of 1 j p (or 2 j p ) indicates that j I is the most likely patch centered at left (or right) NCO.

Fig. 5 .Fig. 6 .
Fig. 5. Visualization of the ONH segmentation produced by our algorithms from 10 randomly selected eyes, compared with A-scan classification-based method and the reference standard.From left to right: NCO detection in central B-scan (number 100 of 200) of the SD-OCT volume (the red dots denote the NCO points, the red lines indicate the disc border and the green line indicate the cup border), A-scan classification-based segmentation, coarse optic disc and cup margin location, patch searching-based segmentation and reference standard.The red and green curves denote the disc and cup margin, respectively.
. The computation cost of an A-scan classification-based segmentation is high because the k-NN classifier cannot use a pre-training model and needs to calculate the distance of each query instance to all training samples.