Active contour method for ILM segmentation in ONH volume scans in retinal OCT

: The optic nerve head (ONH) is aﬀected by many neurodegenerative and autoimmune inﬂammatory conditions. Optical coherence tomography can acquire high-resolution 3D ONH scans. However, the ONH’s complex anatomy and pathology make image segmentation challenging. This paper proposes a robust approach to segment the inner limiting membrane (ILM) in ONH volume scans based on an active contour method of Chan-Vese type, which can work in challenging topological structures. A local intensity ﬁtting energy is added in order to handle very inhomogeneous image intensities. A suitable boundary potential is introduced to avoid structures belonging to outer retinal layers being detected as part of the segmentation. The average intensities in the inner and outer region are then rescaled locally to account for diﬀerent brightness values occurring among the ONH center. The appropriate values for the parameters used in the complex computational model are found using an optimization based on the diﬀerential evolution algorithm. The evaluation of results showed that the proposed framework signiﬁcantly improved segmentation results compared to the commercial solution.

The retina is the only part of the CNS that is readily accessible by optical imaging, putting great potential into its imaging in the context of these disorders. Currently, spectral domain optical coherence tomography (SD-OCT) is the method of choice to acquire retinal 3D images in µm resolution [11,12].
The human retina demonstrates two macroscopic landmarks, whose analysis is especially promising in the context of neurological disorders. The macula around the fovea is the visual field's center and contains the highest density of retinal ganglion cells. Macular SD-OCT images can be analysed quantitatively with intra-retinal segmentation [13]. The derived thickness or volume changes can then be used to quantify neuro-axonal damage in neurological disorders, e.g. in multiple sclerosis [14][15][16].
The second landmark is the optic nerve head (ONH), where all axons from retinal ganglion cells leave the eye towards the brain, thereby forming the optic nerve.
Two membranes define the ONH region and limit the ONH towards the inner and outer eye: the inner limiting membrane (ILM) and the Bruch's membrane (BM) (Fig. 1). Segmenting both structures provides an important starting point for calculating imaging biomarkers of the ONH. Yet, development of suitable segmentation approaches is still an active research topic since ONH segmentation presents several difficult challenges for image analysis. An example where the commercial device fails is shown in Fig. 1. At the ONH, ILM segmentation is particularly challenging because of the broad range of ILM topologies in highly swollen to shallow ONH, ONH with deep cupping, and in some cases even with ILM overhangs. Segmentation is further complicated by a dense vasculature with often loose connective tissue, which can cause ILM surfaces with vastly irregular shapes, see Fig. 10.
A simplified ONH parameter is e.g. the peripapillary retinal nerve fiber layer thickness (pRNFL), which is measured in a circular scan with 12°diameter around the ONH, thereby avoiding direct ONH analysis. Certainly though, the ONH 3D structure itself contains potentially highly diagnostically relevant information as well, as it e.g. has been shown in glaucoma [17][18][19][20]. However, in applications concerning neurological disorders we need not only the minimum distance from the BMO points to the ILM, the so called BMO minimum rim width parameter (BMO-MRW), an already established parameter in the research of glaucoma [17][18][19][20], but also several others that incorporate various regions, volumes and areas. These will be able to capture the changes that ONH undergoes in form of swelling, as for example in idiopathic intracranial hypertension [21,22], as well as in form of atrophy after optic neuritis [23,24] or other optic neuropathies, for example, in the context of multiple sclerosis [25] or neuromyelitis optica spectrum disorders [5].

State of the art and proposed solution
One of the most frequently applied methods for retinal layer and ILM segmentation is based on graph-cuts [13], which works nicely in conditions with smooth surfaces, e.g. at the macula. Despite these important advances, the main drawback of graph-cuts based segmentation is that finding the shortest path through a graph is often implemented via Dijkstra's algorithm, which results are retrieved as a graph function, which can only have one coordinate per axial depth scan (A-scan). This makes it impossible to correctly segment structures passing one A-scan more than once, e.g. by ILM overhangs at the ONH. Another disadvantage is that the employed constraints defined in most of the graph-cuts construction can result in underestimating extreme shapes, which regularly occur at the ONH, e.g. deep cups. An attempt to overcome this problem was presented by [26] by incorporating a convex function in the graph-cut approach, which allowed detection of steep transitions, while keeping the ILM surface smooth. Furthermore, in order to allow the segmentation result to have several coordinates per A-scan, [27] employed a gradient vector flow field in the construction of the columns used in the graph-theoretic approach.
However, these recent advances still do not consider features of the ONH shape like overhangs. As a consequence, a tedious and time consuming manual ILM segmentation correction is necessary before ONH shape parameters can be derived in scans with deep ONH cupping or steep forms [28].
Recently Keller et al. [29] addressed this issue as well by introducing a new length-adaptive graph search metric. This approach is employed though to delineate the retina/vitreous boundary of a different application and region of the retina, namely for patients with full-thickness macular holes.
Novosel et al. [30] proposed a 3D approach to jointly segment retinal layers and lesions in eyes with topology-disrupting retinal diseases that can handle local intensity variations as well as the presence or absence of pathological structures in the retina. This method extends and generalizes the existing framework of loosely coupled level set, previously developed for the simultaneous segmentation of interfaces between retinal layers in macular and peripapillary scans of healthy and glaucoma subjects [31], drusen segmentation in subjects with dry age-related macular degeneration [32] and fluid segmentation in subjects with central serous rethinopathy [33]. Some of these methods are still 2D based, [32,33], and all mainly focus on the macula. Carass et al [34] proposed a 3D level set approach based on multiple-object geometric deformable model that overcomes the voxel precision restrictions that are present in the aforementioned methods. This method also focused on macula scans, and included only images without strong pathological changes.
More recently, deep learning emerged as a possible alternative to established graph-cuts, promising a flexible framework for retina analysis. Fang et al. [35] used a convolutional neural network to predict the central pixel label of a given image patch, and subsequently used the graph based approach to finalize the segmentation, while Roy et al. [36] designed a fully convolutional network (FCN) to segment retinal layers and fluid filled cavities in OCT images, partly based on [37] and [38]. These methods are all based on pixel-wise labelling without using topological relationships between layers or layer shape, which can lead to errors in the segmentation. He et al. [39] proposed a framework to correct topological defects by cascading two FCNs. This method was also applied and tested only on macula scans.
In order to address this issue, we here propose a modified Chan-Vese (CV) based segmentation method with sub-pixel accuracy that is fast, robust and able to correctly detect 3D ILM surfaces regardless of ONH shape complexity or overhangs.
Key features of the developed framework are: • This approach addresses the issue of ILM A-scan overhangs, which are frequently seen in eyes from patients with neurologic or autoimmune neuroinflammatory diseases.
• We include an additional lower boundary constraint in the Chan-Vese method in order to avoid the level set evolution in tissue regions with low contrast.
• The parameters c 1 and c 2 from the level set equation [40] were obtained as a result of an optimization process and further modified after several iteration steps by a scaling factor that incorporates the data locally. This greatly increased the segmentation accuracy.
• We adapted a local fitting intensity energy introduced by [41] into the narrow-band context.
• One of the most important aspects of our work is the thorough investigation and automatic setting of parameters involved in the level set equation through an optimization process. Thus, we optimized and validated our approach using 3D ONH scans from a heterogeneous data set of 40 eyes of healthy controls and patients with autoimmune neuroinflammatory diseases.
• We further validated our segmentation and investigated the influence of image quality on the results using 100 randomly chosen scans from a large database that included 416 ONH volumes of healthy controls and patients with autoimmune neuroinflammatory diseases.
In the following we describe the algorithmic approach as well as its optimization for segmentation performance and computation time in detail.

Region based active contour methods
Active contours have been introduced by [42] as powerful methods for image segmentation. A contour C (also called a snake) is moved over the image domain Ω, where the dynamics are defined by the gradient flow of some suitable energy functional F = F[C]. Here, F[C] depends on the intensity function I = I(x) ∈ [0, 1] of the given gray scale image. Thus, the final segmentation is obtained by finding the energy minimizing contour C for the given image I. In region-based methods, the two regions defined by the contour C are used to model the energy function F, as proposed in a general setting by [43].
Let Ω denote the image domain. In the context of OCT volumes, Ω is 3D, and we are looking for a surface C, that divides the retinal tissue from the vitreous body, i.e. the final contour should be the segmented ILM, satisfying: inside (C) = Ω 1 the retinal tissue (1) outside(C) = Ω 2 the vitreous body (2) The classical CV model [40] approximates the intensity function I(x) by some piecewise constant function with values c 1 and c 2 in Ω 1 and Ω 2 , respectively, x = (x, y, z) is a voxel in the image I. The energy F cv is defined as the weighted sum of a region based global intensity fitting (gif) energy F gif , penalizing deviations of I(x) from the corresponding value c i , i ∈ 1, 2, a surface energy F surf given by the surface area of C, and a volume energy F vol given by the volume of Ω 1 : ds, ds is the area element (5) Note that by minimizing F cv , the surface energy F surf leads to a smooth contour surface C, whereas the volume energy F vol yields a balloon force. Moreover, minimizing F cv with respect to the values of c 1 and c 2 results in choosing c 1 and c 2 as the average intensities in the regions Ω 1 and Ω 2 , respectively. The positive weighting parameters λ 1 , λ 2 , µ, ν control the influence of the corresponding energy term. Finding good values for these parameters is crucial for obtaining the desired results. Usually the parameters λ 1 and λ 2 are taken to be equal and may therefore be set to λ 1 = λ 2 = 1. In section 3.5 it will be explained how values for several parameters of our final model were established.

Challenges in OCT images
The classical CV model [40] is robust in segmenting noisy images, without clearly defined boundaries. Moreover, complicated shapes such as overhangs, various connected components, even topological changes are handled naturally. However, applying the original formulation of CV to OCT scans does not yield good results. As already discussed in the literature, see e.g. [41], CV fails to provide good segmentation if the two delineated regions, in our case, Ω 1 and Ω 2 , are strongly non-uniform in their gray values. Performance gets even worse in the presence of very dark regions inside the tissue, see Figs. 2(a) and 2(b), where a slice (B-scan) of a typical OCT volume scan is depicted, or in regions with extreme high intensities inside the tissue or the vitreous. As a consequence, using local image averages, as proposed in the classical CV model, is not able to provide a satisfactory segmentation.
In the following, we address these problems by adapting the global fitting energy to a local one in the narrow band setting [41] and by using a boundary potential to prevent the contour to move into certain regions. These modifications result in a very stable segmentation algorithm for OCT scans provided that the coefficients weighting the energy terms are chosen suitably.

Global fitting energy
We adapt the global values c 1 and c 2 , obtained after optimization (see subsection 3.5), columnwise (per A-scan) in order to obtain correct segmentation results at the ONH, where the tissue has very low intensity and little contrast to the vitreous, see Fig. 3(a). We scale the values on each column using the formula: where M × N is the size of a B-scan. This significantly improves the segmentation results, see Fig. 3(b). Additionally, in order to prevent the contour to penetrate the retina in regions with dark upper but hyperintense lower layers, see see Fig. 3(c), we rescale these values again, after the interface has almost reached the desired ILM contour. The factor used is computed with the same formula as in Eq. (7), but the maximum intensity considers only voxels from the top of the volume to 77 µm (20 px) below the current interface position. Segmentation results obtained after this scaling step are shown in Fig. 3 We choose to rescale the values of c 1 and c 2 instead of rescaling the column intensities, since the latter would also influence the local fitting energy (lif), presented in section 2.2.3.

Boundary potential
To prevent the contour C, at the ONH, from evolving into dark areas, caused by absence of retinal layers and presence of large vessels, characteristic to this region (see Fig. 2(a)), we introduce a local boundary potential V(x): This potential is set to a very high value ρ at these dark regions, detected as follows: in each column, starting from bottom to top, V(x) is set to ρ until the first gray value I(x) is larger than 45 % of the maximum gray value in that column. All the other voxels are set to 0, see Fig. 2(c) for an example.

Local intensity fitting energy
In images with large intensity inhomogeneities, e.g. caused by varying illumination, it is not sufficient to minimize only a global intensity fitting energy. In order to achieve better segmentation results, [41] introduced two fitting functions f 1 (x) and f 2 (x), which are supposed to locally approximate the intensity I(x) in Ω 1 and Ω 2 , respectively. The local intensity fitting energy functional is defined as: Similar to the global constants c 1 , c 2 , explicit formulas for functions f 1 and f 2 are obtained by energy minimization [41]. Since all calculations are restricted to a narrow band along the contour in our approach, we used compact support kernels for K σ based on a binomial distribution with σ representing the kernel size. The adapted formula for calculating f 1 and f 2 is presented in Appendix 7.2. Note, that these modifications also considerably reduce the computation time as compared to a global convolution.

Our energy model
In a final step, the influence of the local and global intensity fitting energy, are combined similar to the work presented in [41] by introducing another weight parameter ω and to arrive at our modified CV functional F: (10)

Level set formulation
To solve the minimization of the above energy functional, the level set method [40] is used, which replaces the unknown surface C by the level set function φ( Energy minimization is then obtained by solving the following gradient descent flow equation for φ to steady state: F local (x), adapted to our narrow band approach, is given in detail in appendix 7.2. In the above Formula, δ ε denotes the smoothed Dirac delta function having compact support on the narrow band of width ε, see appendix 7.1, and k = div ∇φ | ∇φ | is the mean curvature of the corresponding levelset surface.

Initial contour
For initialization we developed a basic 2D segmentation algorithm to create a start contour (start segmentation). In a first step, morphological filters (erosion and subsequent dilation with 15 × 7 ellipse structure element) and a smoothing filter (Gaussian blur with kernel size 15 × 7 and variance σ x = 6, σ z = 3) are used, to reduce speckle noise. In the next step, we set each pixel with at least 35 % of the maximum column intensity to white. The remaining pixels are set to black. To keep only the tissue connected to the retina, enhanced at the previous step, we set the pixel values of all connected components consisting of less than 2,400 pixels, which corresponds to 0.11 mm 2 in the used OCT data sets, to gray. Finally, in each column the contour is set at the first white pixel from top to bottom. If no white pixel exists, the first gray pixel is taken instead. These processing steps are exemplified in Fig. 4.

Narrow band and reinitialization
Starting with the initial contour, the signed distance function and the reinitialization of the 3D level set function is computed using fast marching [44, 45, s. 86]. We use the L 1 norm to achieve faster calculation time. For the same reason the reinitialization process is done only every 10th iteration step, while the computation for f 1 and f 2 each 5th iteration step. All computation is done on the narrow band around the zero level set of width ε + max(σ x , σ y , σ z ), where σ is the kernel size of K σ , see Eq. (9), ε represents the regularization value for the smooth Heaviside function H ε and Dirac delta function δ ε . Our value for ε is 5.

Removal of non-connected volume parts
A bad initial contour can produce non connected volumes in the tissue and in the vitreous. To improve the segmentation we remove these volumes from the final result, as well as ignore a small area around the image border, where unnatural connectivities might occur due to low contrast or missing image information.

Numerical implementation and computation time
To solve the gradient flow equation, an explicit first order time discretization is used, where for the spatial finite difference discretization we followed [40]. Note, the volume scan has different resolution in the two lateral directions, see below, section 3.5.1. We accounted for this by choosing different spatial mesh sizes h x = h z = 1 and h y = 3. The higher value in y direction reduces the influence of differences between the B-scans, especially near the cup, and reduces computation time. Also, curvature weight parameters, µ x = µ z = 3.5 µ y , are direction dependent in order to create a smooth contour. The lower resolution in y direction is forces the segmentation to include even small tissue that is still connected to the retina. For the chosen parameters see Table 4. The computation time depends on the initial contour and on the shape of the ONH, and takes on average 14.2 s (minimum time is 12.8 s, maximum 15.9 s) on a standard PC (Intel Core i7-5600U (2 × 2.6 GHz) with Debian 9 and gcc 6.3.0). For the ONH scan shown in Fig. 12(a) 15.5 s were needed (0.2 s for creating the initial contour and 15.3 s for segmentation using the proposed CV method). The OCT image size was 384 × 496 × 145 voxels.

Parameter optimization
We used an automatic parameter optimization procedure to find values for the parameters ν, ω, c 1 and c 2 for our computational model (11).

OCT image data
Image data consisted of 3D ONH scans obtained with a spectral-domain OCT (Heidelberg Spectralis SDOCT, Heidelberg Engineering, Germany) using a custom ONH scan protocol with 145 B-scans, focusing the ONH with a scanning angle of 15°× 15°and a resolution of 384 A-scans per B-scan. The spatial resolution in x direction is ≈ 12.6 µm, in axial direction ≈ 3.9 µm and the distance between two B-scans ≈ 33.5 µm. Our database consists of 416 ONH volume scans that capture a wide spectrum of ONH topological changes specific to neuroinflammatory disorders (71 healthy control eyes, 31 eyes affected by idiopathic intracranial hypertension, 60 eyes from neuromyelitis optica spectrum disorders, and 252 eyes of multiple sclerosis patients). We have chosen 140 scans randomly from this database, which presented different characteristics, from scans with good quality up to noisy ones, from healthy but also eyes from patients with different neurological disorders, in order to cover a broad range of shapes. 40 ONH scans were manually segmented and used as ground truth for the optimization process. These will be called the group40. The remaining 100 scans -in the following referred to as the group100 -were used for the validation of the segmentation results and assessment of image quality influence on the segmentation results.
Incomplete volume scans as well as those with retinas damaged by other pathologies were not included.

Error measurement
All 40 scans of the group40 were manually segmented and checked by an experienced grader. From this dataset, 20 images were used for optimization, while the other 20 for validating the results. For one optimization run, 10 files were randomly chosen from the optimization set. The measure used for the minimization process was defined as the sum of the errors for the parameter ν, ω, c 1 and c 2 . An error metric similar to the one described in [46] was employed, where the error is defined as the number of wrongly assigned voxels, i.e. the sum of the number of false positive and false negative. Note that this metric does not depend on the position of the retina. In order to compare different optimization results, the accumulated error of all the 20 scans of the optimization set was used.

Optimization algorithm
The method chosen is the multi-space optimization differential evolution algorithm as provided by GNU Octave [47].
This algorithm creates a starting population of 40 individuals with random values. In our case, an individual represents a parameter set for ν, ω, c 1 and c 2 together with the accumulated segmentation error of the 10 selected volume scans. During optimization, the algorithm crosses and mutates the individuals to create new ones, and drops out newly created or old ones depending on which exhibits larger errors. We allowed for at most 2000 iterations and set box constraints for all four parameters. Note that for each newly created individual, the cost function (error) has to be evaluated by first performing 10 segmentations for the randomly chosen OCT-scans, then calculating the error by comparing with results from manual segmentation. Thus, each iteration step is computationally demanding. The differential evolution algorithm has been chosen since it is derivative free and supports the setting of specific bounds for the parameters. Moreover, we observed a high reproducibility of the finally obtained optimal parameter set. To perform the optimization, we used the Docker Swarm on OpenStack infrastructure from [48], which allowed to do parallel computations on a PC cluster.

Optimization results
In Table 1, the results of all 10 optimization runs sorted in increasing order of the total error are presented. The parameters given in the first line, namely ν = 0.03248, ω = 0.15915, c 1 = 0.24206, c 2 = 0.94945 have been chosen for all subsequent calculations.
Note that the four parameters ν, ω, c 1 , c 2 are not independent from each other as it can be seen in the definition of the energy functional. The parameter that shows the largest variation is the balloon force weight parameter ν. This variation is highly influenced by the presence or absence of one specific volume scan in the randomly chosen optimization set (subset of 10 out of 20), which appears as outlier with highest error in all 10 error distributions as shown in Fig. 5. This occurs because the parameters will account for this particular scan if it is contained in the optimization set.

Evaluation of segmentation performance
Our results are summarized in Fig. 6. For each scan, numbered from 1-40, the first bar (blue) represents the segmentation error given by ILM segmentation implemented in the device, the second bar (green) is the error of our segmentation method starting with an initial contour given by the device's segmentation, and the last bar (red) represents the error of our segmentation method starting with an initial contour calculated as described in subsection 3.1. The lower bars represents the amount of error at the ONH region.
For evaluation and analysis of our results, we used GNU Octave version 4.0.3 [49] and R version 3.3.2 [50]. Our segmentation method outperforms the segmentation available from the device (Wilcoxon Signed-Rank Test p-value = 1.837 × 10 −10 at 0.05 confidence interval).
To check the influence of the initial contour, we compared the results of our segmentation method, obtained starting either with our own initial contour (mean / STD = 21 707 / 5365, STD represents standard deviation) or with the segmentation given by the device (mean / STD = 21 715 / 5152; Wilcoxon Signed-Rank Test p-value = 0.3007 at 0.05 confidence interval). Our model outperforms again the segmentation of the device (error of the segmentation of the device mean / STD = 34 884 / 14 062; Wilcoxon Signed-Rank Test p-value = 1.819 × 10 −12 at 0.05 confidence interval). Moreover, starting with two different initializations, the results of our method are very close to each other (bars 2 and 3), therefore showing that our method is independent on a given rough estimation of the ILM.
We then looked at the error contribution of the central region around the ONH consisting of all A-Scans within a radius of 1.5 mm. Around 50 % of the total error is represented by the errors in the ONH centered region, which usually presents the strongest topological challenges. When evaluating the same comparison with the ILM segmentation computed by the device, our algorithm again performed considerably better (error inside the center region with our algorithm mean / STD = 11 092 / 3306, as compared to the error of the device mean / STD = 18 744 / 6328; Wilcoxon Signed-Rank Test p-value = 1.819 × 10 −12 at 0.05 confidence interval).
Lastly, we evaluated the local error size distribution, i.e. the local Euclidean distance between manual and our automatic segmentation on each B-scan of all 40 scans, see Fig. 7 and mean error as well as 10 %, 50 % and 90 % quantile of the error distribution, see Table 2. Here we have calculated the distance of each line segment of the calculated segmentation line to the ground truth segmentation line at each OCT slice. Note that the manual segmentation is pixel based whereas our proposed method has subpixel accuracy. We can see a large majority of these local  errors are less than 2 µm, followed by a range from 2 µm to 4 µm for most of the remaining ones.

Segmentation validation
To evaluate the performance of our segmentation method in a clinical setting, we analyzed the 100 OCT scans described in section 3.5.1, denoted as group100. An experienced grader manually checked all scans and corrected the segmentation, if necessary. Besides segmentation results, the influence of scan quality on the segmentation was assessed. To this end, another experienced grader labeled regions in all scans (initial group40 see section 3.5.2, and the aforementioned group100) according to the following criteria: • Incomplete scan (top) • Incomplete scan (bottom) • Missing signal • Low SNR • Bad illumination • More than one of the above criteria Figure 8 shows a B-scan with the labeled A-scans that suffer from low SNR, and the corresponding SLO image with a projected map of all A-scans labeled according to these quality criteria. Table 3 shows all A-scans of group100 that have been assessed according to the quality criteria, as well as the number of their corresponding corrected A-scans. It can be seen, that less than 0.2 % of all A-Scans had to be corrected. Note, that scan quality does influence the segmentation results, yet only a smaller percentage of bad quality regions presented also an erroneous ILM: 13.3 % of the A-scans having low SNR and 6.2 % of the ones with low illumination. In the concrete clinical application context a user would be interested on his workload in number of B-scans that must be corrected. The modified A-Scans in Table 3 represent 724 B-scans, a rather  low fraction (5.0 %) of the total number of 14 500 evaluated. Furthermore, Fig. 9 shows the error distribution in µm 3 for all B-scans of group40 and group100, respectively. Although the two groups underwent a different segmentation correction process the results presented in these two figures are similar. Most of the scans present small errors, with few larger outliers. 77.15 % of group40 and 99.3 % of group100 are below 0.0003 mm 3 . To understand the magnitude of these outliers, note that the average volume of a healthy control is (1.80 ± 0.49) [21]. The maximum error per volume for group40, 0.065 46 mm 3 , represents 3.6 % of an average volume, while 90 % of all volumes have an error less or equal to 0.0454 mm 3 , e.g. 2.52 % of the average volume. For group100 the maximum error per volume is 0.012 59 mm 3 , 0.69 % of an average volume, while the error for 90 % of all volumes is less or equal to 0.0025 mm 3 , which represents 0.138 %.
We emphasize again a very important aspect of the ONH OCT scans. Due to the anatomy of the ONH and blood vessels, the contrast between tissue and vitreous is often insufficient to precisely delimit the retinal tissue even for an experienced grader. To investigate the segmentation goodness in these ambiguous parts, a second experienced grader labeled these in all scans of group100. The results showed that, from all the corrected voxels, approximately 25 % were contained in these labeled areas.

ILM segmentation
We present several challenging examples selected from the larger database presented in subsection 4.3. Figure 10 shows a case where a hole in 3D is formed. Our segmentation (red line) is able to correctly detect the ILM, whereas the device implemented method fails. Also Fig. 11 shows several cases where our method manages to detect all visible tissue being connected to the ILM. To emphasize the 3D nature of our method, we present in Fig. 12 rendered surfaces obtained from the segmentation result. It can be clearly seen, that overhangs in the cup create special topologies, that were not properly handled by previous methods.

Discussion
We developed and implemented a method for segmenting the ILM surface from ONH centered 3D OCT volume scans utilizing a modified CV model. We optimized, tested and validated the method on scans presenting a variety of ONH shapes found in healthy persons and patients with autoimmune neuroinflammatory disorders. Our method is providing ILM segmentation that handles not only steep cups, but also challenging topological structures, i.e. overhangs. To achieve this we performed several modifications to the original CV approach, which enabled us to account for low contrast regions especially at the limit between vitreous and retina, for strong inhomogeneities throughout a single B-scan, but also throughout the entire volume by introducing a local fitting energy additionally to the global one. Wang et al. [41] introduced an energy functional with a local fitting term and an auxiliary global intensity fitting one and showed the advantages of using such an approach in terms of accuracy and robustness in inhomogeneous regions. Our method further expands this approach by incorporating this energy formulation in a computational efficient narrow-band method. By locally rescaling the two parameters c 1 and c 2 in the global energy definition (see subsection 2.1 and 2.2.1) we further stop the interface from entering inside retinal tissue. Furthermore, the introduction of a boundary potential successfully accounts for shadows cast by the presence of large blood vessels at the ONH. Not the least, we carefully chose key parameters ω, ν, c 1 and c 2 as a result of an extensive optimization process. As a result, the performance of the proposed multi-energy segmentation framework significantly improved segmentation results compared to the one integrated in the commercial device. Furthermore, our approach's results were independent from the initialization contour.
We did not compare our method with other retinal segmentation methods found in the literature, since each method was evaluated on different types of retinal images, different metrics for the accuracy measurement. However, to our knowledge, none of the existing methods, i.e. when based on graph-cuts, tackles the challenge of segmenting overhangs present especially at the ONH region [27,[51][52][53][54][55][56][57].
Of note, although the length-adaptive method introduced by Keller et al. [29] and evaluated in pathologic eyes in the macula region performed better than the shortest path algorithm, it was not able to successfully segment every feature of the pathology presented. Here, especially weak gradients caused segmentation challenges. Again, a direct comparison with this method is difficult since this approach was applied to the macula region. Also, it is not clear how this approach would handle overhangs that present as several holes in one B-scans.
Several level-sets methods have also been employed in retinal layer segmentation, most of which are related to the macula region and, or specific pathologies in this region [30][31][32][33][34]. [32,33] use 2D segmentation, which limits their capabilities in capturing challenging 3D structures, as seen at the ONH. Furthermore, most of these methods focus on the macula, except [31] which also provides data on ONH scans, although as the authors stated themselves, this algorithm would need further modifications to account for difficult topologies. To overcome the intensity variation and disruptive pathological changes caused by several diseases, [30] proposed a 3D approach steered by local differences in the signal between adjacent retina layers. It is not clear how this approach would perform at the ONH, especially in low contrast between vitreous and interface, since this might be the reason why drusen, especially small ones were not so accurately detected. Another question that arises is, to which extend is the segmentation affected by noise, since the later is highly dependent on the contrast between layers in the attenuation coefficient values.
Deep learning based methods are an important advancement, which might lead to powerful segmentation approaches of retinal tissue [35,36,39]. Their performance is highly dependent on the training dataset [58], and thus their applicability to a broad clinical spectrum of retinal changes remains to be shown. It will be interesting to see first applications of this potentially very promising approach for segmenting ONH scans.
A limitation of our method lies in the successful detection of an initial contour. This failed in one case among the investigated 40 scans used in the optimization process, which was caused by weak image quality and an grossly uneven image contrast. For practical use, a minimum image quality might prevent such errors.
We also tested the performance of our segmentation on a larger database, containing 100 volume scans of healthy controls and several different neurological disorders and investigated the influence of image quality on the goodness of the segmentation. The results have shown that our proposed segmentation is robust against scan image quality (like noise, and bad illumination) with only a small percentage of the total segmentation having to be corrected in these cases. Additionally, overall the amount of error was small, with few outliers, which represented a small amount of the total volume of a complete scan. Therefore, our segmentation represents a considerable improvement over the current commercial solution. Furthermore, it can be employed in the clinical routine.
Our method may present a key component in future approaches to analyze overall ONH structure and volume in health and disease. For this, an accurate detection of tissue between ILM and BM is necessary. Further, detection of the ONH center and the ocular axis are required to establish coordinate origin and scan rotation. We previously developed a 2D based segmentation approach for ONH volume computation [59], which was able to robustly detect the BM as lower boundary in healthy people as well as in patients suffering from various neurological disorders [21,22,60]. However, in this method the other key surface for volume computationthe ILM -was taken from the device itself, which led to erroneous boundary definitions in many cases. Thus, with the proposed modified CV based method, ILM detection for ONH volume computation may be further optimized.

Conclusion
It is crucial to have a robust ILM segmentation for the analysis of ONH structure in both neurologic and ophthalmologic conditions. A better insight in morphometric changes of the ONH could potentially improve diagnosis and the understanding of diseases like multiple sclerosis, neuromyelitis optica spectrum disorders, Parkinson's disease and Alzheimer's dementia.
Our approach provides a flexible and accurate computational method to segment the ILM in 3D OCT volume scans. We could show the superiority in performance of the proposed CV model compared to a standard commercial software method and its robustness against the large variations in the ONH topology. Furthermore, our method is fast and allows for flexible initialization. Thus, it shows great potential for clinical use as it is capable to robustly handle ONH data of both healthy controls and patients suffering from different neurological disorders. We hope that similar concepts, especially the optimization framework, will be applied in other applications as well.
We made this choice in order to have compact support of δ ε , being the narrow band of width 2ε.

Convolution with compact support kernel
The calculation of the local intensity fitting force involves convolutions which have been adapted for allowing to consequently use the narrow band. Here, not all values are well defined and therefore the convolution formula has to be modified. Now the local fitting forces are calculated as follows: If ∫ f 1 (y) nan K σ (y − x) > 0 and ∫ f 2 (y) nan K σ (y − x) > 0, then in other case, set F local (x) = 0 (18)

Parameters
All parameter values of our computational model are shown in Table 4.