1. 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, age-related 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–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–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].

 

Fig. 1 One B-scan of ONH-centered SD-OCT image. (a) The normal case; (b) Glaucoma case.

Download Full Size | PDF

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].

 

Fig. 2 ONH centered SD-OCT B-scan from a normal eye with 9 manually segmented surfaces defining 8 retinal layers.

Download Full Size | PDF

The contributions of this paper are summarized as follows:

2. 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.

 

Fig. 3 Flowchart of the proposed method.

Download Full Size | PDF

2.1 Shared-hole graph search with locally adaptive constraints

2.1.1 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–39]. In the graph search algorithm, volumetric image is defined as a 3D matrix $I(x,y,z)$ with size $X\times Y\times Z$, and the surface is defined by a function $z=f(x,y)$, where $x\in \left\{0,\cdot \cdot \cdot ,X-1\right\}$, $y\in \left\{0,\cdot \cdot \cdot ,Y-1\right\}$, and $z\in \left\{0,\cdot \cdot \cdot ,Z-1\right\}$.

For single surface segmentation, two smoothness constraints, ${\Delta }_x$ and ${\Delta }_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 node-weighted 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

 

Fig. 4 Graph construction for graph search with surface smoothness constraints (a) and surface separation constraints (b).

Download Full Size | PDF

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. $\text{K}$ sub-graph $G_i(V_i,E_i)(i=1,\mathrm{....},k)$ are constructed for$k(k\ge 2)$ surfaces segmentation. Inter-surface arcs are enforced to restrict the maximum and minimum horizontal distances between surface pairs. Figure 4(b) illustrates the surface separation constraints for a pair of surfaces. The maximum surface separation constraint is set as ${\delta }_{i,j}^{\max }(x,y)=2$, and minimum surface separation constraint is set as ${\delta }_{i,j}^{\min }(x,y)=1$.

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.

2.1.2 Shared-hole graph search

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.

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 intra-column 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 shared-hole 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.

 

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.

Download Full Size | PDF

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.

2.2 Preprocessing

2.2.1 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.

2.2.2 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,

Figure 6 shows an example of the results obtained by the proposed method, Otsu thresholding, and conventional graph search method with small constant smoothness constraints (${\Delta }_x=4$ and ${\Delta }_y=25$) or big constant smoothness constraints (${\Delta }_x=25$ and ${\Delta }_y=25$). 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 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.

Download Full Size | PDF

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. Figure 7(b) shows the result after B-scan alignment.

 

Fig. 7 The y-z image before (a) and after (b) B-scans alignment. The highlight box region indicates the distortion region.

Download Full Size | PDF

2.2.3 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 B-scans, 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.

 

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.

Download Full Size | PDF

2.3 Optic disc boundary detection

The optic disc boundary detection method is described in this section. The method begins from several intra-retinal surface segmentations, which is used for feature extraction. Then a random forests classifier is trained to classify the optic disc region. After classification, the result is further refined to reduce the interference of vessel shadows.

2.3.1 Intra-retinal surface segmentation

In this step, four surfaces, namely surface 1, 6, 8, and 9 are segmented with the proposed graph search method with locally adaptive constraints. Note that surface 1 is already detected in original B-scans, and its position can be polar-transformed to the corresponding position in radial B-scans. Then the radial B-scans are down-sampled by a factor of 2 both in z-direction and x-direction to obtain a low-resolution image volume. Surface 9 is first detected in this volume and then refined in the original volume. The locally adaptive constraints for surface 9 detection are learned from the polar-transformed surface 1 (as $f_{S1}()$) in the same resolution, similarly as Eqs. (4) and (5), only replacing $f_{Otsu}()$ with $f_{S1}()$. The default smoothness constraints for graph search are set to $d_x=1$, $d_y=1$ empirically, and the search regions are limited to the region below detected surface 1. For the refinement, graph search is confined in the volume centered at the interpolation of the initial low-resolution result and with height of 20 voxels.

Similarly, for surface 6 and 8, a double surface graph search is applied with constraints learned from surface 9 in low-resolution volume and refined in the original volume.

2.3.2 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,

Table 1. Features for Optic Disc Region Classification

View Table | View all tables in this article

 

Fig. 9 41 feature maps extracted from a test data.

Download Full Size | PDF

 

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).

Download Full Size | PDF

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.

 

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

Download Full Size | PDF

2.4 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.

 

Fig. 12 The down-sampling process.

Download Full Size | PDF

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 $d_x=1$, $d_y=1$ empirically.

3. Experimental method

3.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\times Y\times Z$).

3.2 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.

3.3 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:

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 B-scan, 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.

4. Results

4.1 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.

 

Fig. 13 Permutation accuracy importance for 41 features.

Download Full Size | PDF

Table 2. The DSC value of optic disc boundary detection (Mean ± SD)

View Table | View all tables in this article

 

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 B-scan 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].

Download Full Size | PDF

4.2 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 3. Average Unsigned Border Positioning Error of ILM Surface Segmentation

View Table | View all tables in this article

 

Fig. 15 Example result of 9 surfaces segmentation on the B-scan next to the ONH region and on the center B-scan of ONH region. (a), (b) Segmented surfaces by conventional multi-resolution graph search with constant constraints [44]. (c), (d) Segmented surfaces by Zang’s method [30]. (e), (f) Segmented surfaces by OCTExplorer software [53]. (g), (h) Segmented surfaces by GS-AAM method [31]. (i), (j) Segmented surfaces by proposed method. (k), (l) Ground truth.

Download Full Size | PDF

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.

Table 4. Summary of Mean Unsigned Border Position Error for All Data

View Table | View all tables in this article

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.

4.3 Computational time

The proposed algorithm was implemented in C++ and tested on a PC with Intel i7-6700 CPU@3.40GHz and 32GB of RAM, where only single core was utilized. The average running time of the algorithm for optic disc boundary detection and 9 layers segmentation in one 3D volume data was 263.9 seconds.

5. Discussion

In this paper, we have proposed a novel method for automated segmentation of retinal layers and optic disc region on ONH centered SD-OCT images both for glaucomatous eyes and normal eyes. The optic disc region was first detected by random forest method on polar-transformed radial OCT B-scans with features representing both textural and structural information. A mean DSC value of 0.925 was achieved, which was statistically significantly higher than Hu’s method [25] and Zang’s method [30] both with p-value <0.001. We also trained 4 random forest models with feature 1-9, feature 10-41, top 30 and top 35 features by ranked permutation accuracy importance, respectively. And the DSC of the optic disc boundary detection results obtained from each model were compared (Table 2), which indicated the effectiveness of the 41 features. The existence of border tissues confused Hu’s method and caused erroneous optic disc boundary detection, and Zang’s method [30] failed to detect optic disc boundary accurately, due to the strong noise in the test data and the changeable ONH structure among glaucoma and normal eyes, but the proposed classification method with novel features can handle these conditions better.

Then two layer segmentation methods, namely Otsu segmentation guided graph search method and shared-hole graph search method, both with locally adaptive constraints were proposed for 9 surfaces segmentation. The former method was used for segmenting internal limiting membrane surface. Due to deep cupping shape inside the ONH region, the multi-resolution graph search method [44] with global constant surface smoothness constraints, Zang’s method [30], GS-AAM method [31] and OCTExplorer software [53] all generate erroneous segmentation result in most cases. To solve this issue, the Otsu segmentation guided graph search method learned the surface smoothness constraints from Otsu thresholding result, confined the searching space near the thresholding-based result, and finally achieved segmentation by minimizing the gradient-based cost function. Thus, the method combined the merits of both region-based and boundary-based segmentation methods. By comparison, the proposed method had the best segmentation result with mean unsigned border errors of 5.38 ± 4.23 in μm.

The latter method was designed for segmentation of the rest 8 surfaces, for which the existence of hole structure (optic disc region) hindered the segmentation accuracy. The proposed shared-hole graph search method effectively eliminated the effect of optic disc region by adding specific arcs in the graph construction step. Locally adaptive constraints were learned from previously segmented layers, and used as constraints for the surface smoothness for the subsequent layer segmentation. Compared with other 4 reference methods, our proposed method performed better with the mean unsigned border errors of 7.51 ± 5.49 in μm for the rest 8 surfaces.

In most cases, the speckle noise in OCT images can strongly hamper optic disc detection and layer segmentation. So, effective denoising methods such as [50–52] are needed to improve the quantitative analysis performance.

The idea of locally adaptive constraints for graph search can also be applied in segmentation of macular centered OCT scans both for normal and abnormal cases. As the smoothness constraints are learned from the data itself, off-line training from other data sets is not necessary, and bias from the training data is avoided. In the future, we will apply the method for layers segmentation in OCT scans of eyes with pathologies such as choroidal neovascularization (CNV) and drusen.

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.

6. 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.