Optic nerve head three-dimensional shape analysis

Abstract. We present a method for optic nerve head (ONH) 3-D shape analysis from retinal optical coherence tomography (OCT). The possibility to noninvasively acquire in vivo high-resolution 3-D volumes of the ONH using spectral domain OCT drives the need to develop tools that quantify the shape of this structure and extract information for clinical applications. The presented method automatically generates a 3-D ONH model and then allows the computation of several 3-D parameters describing the ONH. The method starts with a high-resolution OCT volume scan as input. From this scan, the model-defining inner limiting membrane (ILM) as inner surface and the retinal pigment epithelium as outer surface are segmented, and the Bruch’s membrane opening (BMO) as the model origin is detected. Based on the generated ONH model by triangulated 3-D surface reconstruction, different parameters (areas, volumes, annular surface ring, minimum distances) of different ONH regions can then be computed. Additionally, the bending energy (roughness) in the BMO region on the ILM surface and 3-D BMO-MRW surface area are computed. We show that our method is reliable and robust across a large variety of ONH topologies (specific to this structure) and present a first clinical application.


Introduction
Many neurological and ophthalmological conditions affect the retina. Two landmark structures of this affection are the macula and the optic nerve head (ONH). Parameters describing the shape of these structures, especially 3-D parameters, are gaining more and more importance in understanding normative data and disease-specific changes, in form of additional information along with traditional quantitative ones. For example, in the case of the macula, several studies and approaches have already been published in order to provide a better insight into the definition of healthy structure and shape. [1][2][3][4][5] ONH represents the part in which all retinal ganglion cell axons gather to leave the eye toward the brain. As such, investigating the ONH is central to the diagnosis of many disorders affecting the optic nerve like glaucoma, 6 idiopathic intracranial hypertension (IIH), 7,8 optic neuritis (ON), 9,10 and optic neuropathies of other etiology, for example, in the context of multiple sclerosis (MS) or neuromyelitis optica spectrum disorders (NMOSD). [11][12][13] The advent of fast spectral domain optical coherence tomography (OCT) recently allows in vivo 3-D imaging of retinal structures, including the ONH. While high-resolution ONH images can now be taken in 3-D, the complex ONH shape has so far made automatic image analysis challenging.
Main focus of this work is to fill in this gap by presenting and validating a robust and fully automatic method that is capable of analyzing ONH shape over a range of conditions. Our approach includes a fully automatic 3-D shape analysis pipeline, with various techniques that will be presented throughout the paper. The algorithm is shown to be stable to various 3-D ONH scans and is-to our knowledge-the first approach that fully exploits the 3-D nature of the ONH for computational image analysis.

Related Work
Analyzing tissue damage and structural changes in the ONH is one of the key goals to improve the diagnosis and understanding of diseases related to this structure. The main focus of ONH research lies in the field of ophthalmology, the most prominent topic being glaucoma. The most common parameters utilized include length, areas, volumes, or ratios to quantify various regions of the ONH. The main anatomical structures needed to compute these parameters are the inner limiting membrane (ILM), Bruch's membrane or the lower boundary of the retinal pigment epithelium (denoted throughout the paper, for simplicity, by RPE) and the Bruch's membrane opening (BMO) points. The first two structures comprise the retina, and the BMO points have recently gained more and more attention since Reis et al. 14 introduced these as being the true anatomical structure defining the optic disc as to the clinically identified margin using fundus photography. Following this work, several other groups 15,16 proved BMO-based parameters superiority in reliability compared to cup-to-disk ratio.
For ONH analysis, Enders et al. 17 and Muth and Hirnei 18 used manual segmentation of the BMO points and minimal rim width (BMO-MRW), as well as commercially available ILM surface detection with manual correction to study the correlation between visual fields and structural ONH changes. Chauhan et al. 19 proposed BMO-MRW as a marker for early glaucoma detection also using manual segmented BMO points and BMO-MRW measurements, while Pollet-Villard 20 used a similar manual process to prove that structure-function relationship was significantly stronger using BMO-MRW over other ONH parameters derived from spectral-domain OCT. Furthermore, several studies focused on the description of normative values that characterize the shape of the ONH. Chauhan et al. 21 analyzed BMO-MRW, the orientation of the long axis of BMO in a normal population, while Enders et al. 17 described BMO-MRW in micro-and macrodiscs. 22 BMO minimum rim area (BMO-MRA) correlates with the total number of nerve fibers traversing the optic nerve and was introduced and investigated in Refs. 23 and 24. Both studies have shown BMO-MRA's high diagnostic power for glaucoma. All these studies used commercial available OCT software to retrieve the ILM surface and BMO points, which were then manually controlled and corrected by experienced graders.
Another important parameter is the ONH volume, previously analyzed mainly in relation to conditions characterized by ONH swelling. To this end, a 2-D segmentation has been previously developed by our group. The method 25 is able to robustly detect the lower boundary of the RPE in healthy as well as in ONHs from patients suffering from various neurological disorders that lead to swelling of the ONH. Several other publications investigating the same condition 26,27 followed, all using a graph-cut based approach to segment the ILM and BM at the ONH. Lee et al. 28 introduced an end-to-end pipeline to compute several morphometric parameters (volumes in different ONH regions and related area parameters). Furthermore, a surface correspondences approach to create a normalized space was presented. However, this pipeline is semiautomatic and needs manual segmentation of the ILM, BM, and BMO points. For the development of the mean surfaces, an accurate registration between the ONH surfaces is needed. Gibson et al. 29 introduced an ONH registration algorithm to compute the one-to-one correspondence between two ONH surfaces using the hemispherical surface and volume registration. Later, Lee et al. 30 proposed a more sophisticated registration algorithm based on surface currents and hemispherical demons. Recently, the shape variability of retinal nerve fiber layer (RNFL)-choroidal thickness was investigated using a nonrigid surface registration for longitudinal analysis. 31 In general, most of the published methods related to ONH morphometry are computing the traditional parameters (lengths, widths, and volumes) and performing layers segmentation either manually or semiautomatically. To the best of our knowledge, none of the published ONH morphometry methods use proper manifold 3-D ILM or BM surfaces but rather a 2.5D surface (i.e., a graph function on an XY-grid).

Contribution
In this paper, we propose a fully automatic 3-D shape analysis of the ONH region and introduce 3-D shape parameters along with commonly used ones. The objective of our proposed method is to retrieve robust and reliable 3-D quantitative measurements that describe different aspects of the various shapes of the ONH. Specifically: • We introduce a fully automatic pipeline for 3-D shape analysis of the ONH region.
• We automatically segment ILM and the lower boundary of RPE along with BMO points.
• We compute triangulated 3-D surfaces of the ILM and the lower boundary of RPE.
• We introduce several 3-D shape analysis parameters, along with traditional parameters.
• We prove the robustness of our algorithm for repeated measurements in healthy control (HC) data.
• We provide quantitative parameters for different ONH regions in healthy people with normal vision, as well as in a dataset with clinical relevance for ONH swelling and atrophy.

Method
In this section, we explain the procedure to compute several 3-D shape parameters of the ONH including the preprocessing of the ILM and the RPE surfaces and correspondence between them. Figure 1 shows the algorithm pipeline, where an OCT volume scan is the input of the algorithm. Next, we compute the triangulated ILM and RPE surfaces along with the BMO points. These three structures represent the ONH 3-D shape and serve as inputs for further shape analysis of the ONH. In the remainder of the paper, the ILM and the RPE surfaces will be represented as M ILM and M RPE , respectively, and the BMO points are denoted by P. The ILM and the RPE surfaces are triangulated manifold surfaces and can be written in terms of the set of vertices and faces (triangles): E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 ; 3 2 6 ; 3 5 0 M ILM ¼ fV ILM ; F ILM g and M RPE ¼ fV RPE ; F RPE g: These two surfaces can have different numbers of vertices and triangles. Let us consider that n ILM and m ILM represent the numbers of vertices and faces in the ILM surfaces. Similarly, n RPE and m RPE denote the size of V RPE and F RPE , respectively. The BMO points are represented as P ¼ fp i ∈ R 3 ji ¼ 1 · · · n p g, where n p is the number of the BMO points.

OCT Image Data
Our algorithm has as input 3-D ONH scans obtained with a spectral-domain OCT (Heidelberg Spectralis SDOCT, Heidelberg Engineering, Germany) using a custom protocol with 145 B-scans, focusing the ONH with a scanning angle of 15 deg ×15 deg and a resolution of 384 A-scans per B-scan. The spatial resolution in x direction is ≈12.6 μm, in Fig. 1 The pipeline of the proposed algorithm. axial direction ≈3.9 μm, and the distance between two B-scans is ≈33.5 μm.

Smoothing and intensity normalization
The RPE lower boundary represents the termination of the retina and is therefore an important parameter in several morphometric computations. In this subsection, we present our RPE segmentation approach based on the method presented in Ref. 25. First, several preprocessing steps are performed, often employed in OCT images. Consider Iðq xy Þ the intensity of a pixel q xy . Our algorithm starts by applying a large Gaussian smoothing filter (σ ¼ 5 pixels isotropic with kernel size ¼ ð10 μm × 14 μmÞ on each B-scan separately. The smoothing operation not only reduces speckle noise present in most OCT data but also facilitates the approximation of the two most hyperintense layers, the RNFL and the RPE. Then, to address varying intensities, a contrast rescaling on each slice (B-scan) is performed. Contrast inhomogeneities can occur in the form of a B-scan having regions with different illumination or as several B-scans of the same volume with very different intensity ranges. Specifically, a histogram-based amplitude normalization method 32 is used to map the signals in the original image linearly between [0, 1] using as low cutoffs the first, 66th percentiles and as high cutoff the 99th percentile on the histogram of the B-scan, where the sampled column (A-scan) is located. Figure 2(b) shows one B-scans of the resulting smoothed and normalized original volume. Figure 2(a) shows the same B-scan with its original gray values.

RPE approximation
First, we start by approximating ILM as the upper boundary. At each A-scan, we detect the first pixel from top in the smoothed and normalized volume, I SN , higher than 1/3 of the maximum value in the B-scan containing the A-scan. This gives us a set of initial estimate points for the ILM, denoted by ILM init . Next, we approximate the upper boundary of RPE. First, we compute the image derivative, ∇I SN , of each B-scan (vertical gradient) using a Sobel kernel. Looking along each A-scan, starting from the ILM init set, we perform several intermediate steps for the RPE upper boundary approximation. We estimate inner and outer segment junction regions (ISOS) by finding the first set of points p, as shown in Fig. 2

(c):
E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 3 2 6 ; 6 0 2 p xy ¼ f∇I SN ðq ixy Þ < ð∇I SN ðs xy Þ∕3Þg; s xy ∈ ILM init : (1) Starting from the set p xy , RPE upper boundary is approximated: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 5 4 9 RPE upper ¼ max½∇I SN ðp xy Þ; (2) where we only consider the points below IS/OS: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 5 0 6p xy ∈ fks ixy − p ixy k < 60 μmg: At this stage, our input is a list of points that belong to the upper boundary of RPE in each B-scan. This list comprises among the points correctly positioned at the upper RPE boundary, also several outliers, especially in the presence of shadows cast by blood vessels, as well as at the region of the ONH. In order to remove these, we make use of the observation that the first and last third of each B-scan most probably lie outside the optic nerve hear area. This is an observation valid for all the ONH centered scans independently of device or scan settings. Even in these two parts, outliers caused by shadows of blood vessels might be present. To remove these, we take the gradient of the line consisting of the position of upper RPE points and compute the mean value of these from coordinates that most likely belong to the correct upper RPE points. These coordinates represent RPE boundary points that form the largest part of the gradient line between outliers (outliers in the gradient are considered to be >40 μm). The first seed point is then detected as the one closest to the mean value. Starting from this seed q seed , we iteratively remove outliers from RPE upper (points where kq seed − q new k > 70 μm). Analog outliers from the last third of the B-scan are removed. The resulting point set of one Bscan is shown by the blue points in Fig. 2 Correction of these errors is done by applying a cubic smoothing spline with a high smoothing parameter. Note that in the case of presence of blood vessels, large regions of missing coordinates for the RPE might occur and the cubic spline can present deviations from the desired smooth contour.

TPS fitting to the RPE lower boundary surface
Finally, to account for motion artifacts in consecutive B-scans, but also for the natural curvature of the retina, we propose an efficient two-stage thin-plate spline fitting (TPS), which improves the approach proposed by Ref. 33 without making use of the orthogonal scans presented in the work of Ref. 34. First, TPS least-square approximation is performed. The number of control points used is determined by the size of the surface along each axial dimension. At this stage, the number is set to 25% in the slow scan direction, 15%, respectively, and the control points are evenly distributed along each direction. This enables the TPS, in combination with a smoothing parameter, α, see Ref. 35 set to 0.85, to create a more smoothed surface, which keeps the curvature of the retina without being influenced from motion artifacts especially along the slow axis. In our experiments, we found that values of α ¼ ½0.70; 0.85 provide consistent results. Extreme grid points in the original surface defined as mean + standard deviation in local nonoverlapping neighborhoods of 10 × 10 grid points of the TPS surface are removed. Then, the actual TPS fitting similar to Ref. 34 is applied. The choice of parameters at this second step is strongly influenced by the fact that we reduced outliers at a previous stage. Specifically, for grid points, we use 20% in the slow scan direction, 10% in the fast scan direction, with smoothing parameter 0.45. Consistent results were obtained for α ¼ ½0.30; 0.50. Our strategy is to obtain a TPS closer to the data in the grid points while smoothing the artifacts present in the position of the detected RPE lower points, especially at the presence of blood vessels, or in the close vicinity of the approximated ONH region. Both stages are done on the RPE lower without including A ONH . Figure 3(a) shows the original RPE surface with typical artifacts in the in-between B-scans direction. These are corrected after applying our TPS approach while keeping the shape of the original surface. The result is presented in Fig. 3(b). The result of the RPE lower boundary after performing the TPS is shown in Fig. 2(e).

BMO Points Computation
BMO is the termination of the Bruch's membrane (BM) layer and was proposed as a stable zero reference plane for ONH quantification. 14 It is a key parameter in the detection of ONH shape parameters. A challenge in BMO detection is the correct identification of these points, especially in the presence of shadows caused by blood vessels, or the border tissue of Elsching, 14 a structure similar to the BM. We propose a BMO points'

Volume flattening
We start by flattening the whole OCT volume. This step, performed in several retinal segmentation algorithms, 33,36 refers to the translation of all A-scans such that a chosen boundary in the volume is flat. We choose to align the retina to the smoothed RPE lower . The alignment facilitates the volume reduction process, as well as the differentiation of BMO from other tissue.

Volume reduction and vessel suppression
We reduce the OCT volume to a region comprising only the BMO. Depending on the ONH region fitted by the ellipse and the position of RPE lower , the reduced volume can change from 384 × 496 × 145 voxels in its original size to 384 × 90 × 52 voxels, by taking only the B-scans containing the ellipse and in z direction RPE lower AE 100 μm. Vessels appear in the RPE layer as dark intensity regions or shadows. These affect the detection of the BMO, especially because vessels gather at the ONH creating large shadows. Our approach focuses on emphasizing the RPE layer and suppressing these artifacts.
To this end, we first smooth the reduced OCT volume containing the original gray values with an anisotropic diffusion filter 38 and then apply a 2-D Morlet wavelet filter 39 for each B-scan to enhance the RPE line.

BMO points detection
The end-points of the rough ONH area, A ONH , provide the starting points for BMO points detection. On each B-scan, the starting points are updated with new BMO points candidates if they meet the following conditions: (1) have minimum value in the 2-D Morlet filtered image, (2) dðp new ; p BMO-seed Þ < 30 μm, and (3) the curvature in a neighboring region of five voxels is almost 0 to avoid including the tissue of Elsching. In case the BMO points detected in the left and right part of one B-scan overlap, the BMO starting or end region previously defined by A ONH are updated accordingly. An example of a pair of (left and right) detected BMO points are shown in Fig. 2(f).

ILM Surface Computation
The ILM separates the retina from the vitreous body and defines a critical boundary layer for the ONH. Several ILM segmentation methods have been published to separate the ILM layer around the ONH region, and most of them compute the ILM layer as a graph function and are unable to capture complex and variable topological structures of the ONH. To compute the ILM surface, we previously developed a method introduced by Gawlik et al. 40 This method is based on an active contour method of Chan-Vese type and produces a truly 3-D ILM segmentation unlike other state-of-the-art methods. Figure 4(a) shows the ILM surface, which is computed using the marching cubes algorithm, 41 where the input level sets are computed using the method proposed by Gawlik et al. 40 The lower right corner of Fig. 4(a) shows that the approach is capable to reconstruct a complex topological structure in the ONH region and the marching cubes algorithm produces a manifold ILM surface M ILM with proper neighborhood and properly oriented face normals.

ILM Surface Smoothing
During data acquisition of the ONH region using the 3-D OCT scanner, noise is inevitable due to various internal and external factors. As it can be seen from Fig. 4(a), the ILM surface is not smooth and has various noise components. Staircase artifacts are also shown in the left corner of Fig. 4 as an effect of the marching cubes algorithm and of the volume scan resolution. For an accurate shape parameter computation, these artifacts including noise components should be removed at the early stage of the algorithm. To compute a noise free and a high-fidelity ILM surface, we use a robust mesh denoising algorithm, proposed by Yadav et al. 42 In general, mesh denoising algorithms are divided into two categories: isotropic and anisotropic methods. Isotropic methods remove noise effectively but produce a volume shrinkage, which leads to an incorrect shape analysis of the ONH. Anisotropic methods are feature preserving mesh denoising algorithms and induce a small volume shrinkage compared to the isotropic methods. In the case of the ILM surface, the anisotropic methods treat the marching cube artifacts as features and lead to an incorrect shape analysis. The method in Ref. 42 is a combination of both isotropic and anisotropic methods and produces a high fidelity smooth ILM surface without staircase artifacts and minimum volume shrinkage, as shown in Fig. 4(b). This method also produces a high quality triangular mesh with proper face normal orientation, which is vital in further shape analysis of the ONH.

Ellipse Fitting to BMO Points
As presented in Sec. 2.3, the proposed algorithm computes the BMO points automatically. Due to blood vessels around the ONH, noise components and 3-D OCT scan patterns, the BMO points are nonuniform and noisy, as shown in Fig. 5(a).
To remove these artifacts, we fit an ellipse to the BMO point onto XY-plane. First of all, the BMO points P are projected onto the corresponding XY-plane and denoted as a points set: P 2D ¼ fp i ∈ R 2 ji ¼ 0; · · · ; n p − 1g. Then, we apply the Another key parameter in the ONH shape analysis is the center of the BMO points. This is computed as the barycentric of the all P 2D points: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 6 3 ; 5 2 1p wherep c ∈ R 2 . Figure 5 shows that the ellipse fitting is not only removes the noise but also increases the data points uniformly.

Correspondence between ILM and RPE Surfaces
The total retina at the ONH region is delimited by ILM and RPE. For further analysis, it is necessary to find the corresponding points between these two surfaces. In the proposed method, we compute vertices in the RPE surface corresponding to each face (f i ∈ F ILM ji ¼ 0; · · · ; m ILM − 1) of the ILM surface.
In general, the RPE surface, represented here as a function graph: M RPE ∶R 2 → R, has a less complex structure compared to ILM. In the OCT scanner, the number of the A-scans (x-direction) and the number of the B-scans (y-direction) are fixed, which create a regular XY-grid as a domain for the RPE graph function. Therefore, the index of each vertex of the RPE surface can be computed using the number of x-lines (vertical lines) and y-lines (horizontal lines) and the sampling size in both directions, denoted by ε x and ε y , respectively. The numbers of x-lines and y-lines are computed using the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 6 3 ; 2 3 3 where d e represents the ceil function and i denotes the corresponding vertex in the RPE surface. For an accurate computation, we check the neighborhood of vertex i of the RPE surface and the corresponding vertex is computed as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 7 ; 3 2 6 ; 6 2 1ṽ i ¼ṽ j ∈ Ω i j min jc i −ṽ j j; where Ω i represents the 3 × 3 neighborhood (at XY-plane) of vertex i. The termṽ i is the projection of vertex v i ∈ V RPE onto XY-plane. Finally, we get the set C ¼ fv i ∈ R 3 ji ¼ 0; · · · ; m ILM − 1g, which represents the set of RPE surface vertices corresponding to each face in F ILM . A visual representation of the correspondence computation is shown in Fig. 6, where the RPE's regular XY-grid is shown in blue and the projected ILM's vertices and edges are painted in red.

BMO Region Segmentation
For the ONH shape analysis, the region inside the BMO points is of special interest since BMO points represent the optic disc margin. To segment this region, we exploit the elliptic representation of the BMO points in R 2 along with their center, as mentioned in Eq. ilmc ¼ j ∈ 0; · · · ; n ILM − 1j min jṽ ILM j −p c j; whereṽ ILM j ∈Ṽ ILM andṽ rpe j ∈Ṽ rpe . The vertices v ILMc ∈ V ILM and v RPEc ∈ V RPE represent the center of the ILM and the RPE surfaces, respectively.
To compute the BMO regions in the ILM and the RPE surfaces, we transform the ellipse into a circle using the affine transform. Let us considerp e i represents the fitted ellipse point top i as mentioned in Sec. 2.6 andp c i denotes the corresponding affine transformed point on a circle of radius r. This transformation reduces the complexity of the BMO region computation and  Journal of Biomedical Optics 106004-6 October 2018 • Vol. 23 (10) improves the speed of the algorithm. Now, by using the circular representation of the BMO points, the BMO regions in both RPE and ILM surfaces are computed as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 6 3 ; 7 1 9 where r is the radius of the transformed circle andc ILM j ,c RPE j represent the centroid of a face in the ILM and the RPE surfaces, respectively. The computation shown in Eq. (9) is done using the disk growing method.

ONH cup volume
The ONH cup is defined as a segment of the ILM surface, which is below the RPE surface, as shown in Figs. 7(a) and 7(b). Note that the cup is not present in every ONH volume scan. For example, mostly in the case swollen ONH scans, the ILM surface is always above the RPE. To detect the availability of the ONH cup, we go to each face f i ∈ F ILM and compute its centroid c i . As we discussed in Sec. 2.7, each face f i of the ILM has the corresponding vertex v i ∈ C in RPE. If c z i − v z i ≥ 0 for all faces in ILM, then, there is no cup available in the ONH region. Otherwise, there is a cup. The terms c z i and v z i show the corresponding z-coordinates (height). Similarly, we can compute the cup region: where Ω cup consists of faces (triangles) of ILM, which are below the RPE surface. As it can be seen from Fig. 7(b), the cup region is also a manifold surface with proper face normal orientation.
To compute the volume of the cup accurately, we exploit the face normal information at each triangle of the region. The cup volume is computed using the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 6 3 ; 3 4 3 V cup ¼ where A i represents the area of a triangle, which is a projection of the face f i of the ILM surface onto XY-plane and h i is the height with respect to the RPE surface. These variables are defined as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 2 ; 3 2 6 ; 7 5 2 whereẽ 0 andẽ 1 are the connected edges of the projected triangle. The cross-product between the two edges will take care of the orientation of the corresponding face and enables a precise volume computation even in complex topological regions.

Central ONH thickness
The CONHT is defined as the height difference between the center of the ILM and the RPE surfaces, as shown in Fig. 7(c). The CONHT of ONH volume scan is computed as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 3 2 6 ; 5 9 0 where v z ilmc and v z rpec show the corresponding height value (z-coordinates) of the center vertices v ilmc and v rpec , respectively, introduced in Sec. 2.8.

BMO region volume
The BMO region volume is computed using the segmented ILM and RPE surfaces. We separated the cup volume from the BMO region volume such that it does not include the cup volume, if it exists. Then, the BMO region volume can be computed as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 3 2 6 ; 4 3 9 Similar to the ONH cup volume, we compute the BMO region volume using the similar formula: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 3 2 6 ; 3 8 5 where A i is the area of the face f i that belongs to the set Ω BMO and is computed using Eq. (12). Similarly, h i is also computed using the corresponding vertices in the RPE surface as mentioned in Eq. (12).

ONH total volume
Similar to the ONH cup and the BMO region volumes, the ONH total volume is also computed using ILM and RPE surfaces. The total volume is computed from the circular region, with radius 1.5 mm, centered at v ilmc and v rpec for ILM and RPE surfaces, respectively, as shown in Fig. 8(a). Similar to Eq. (9), we compute the circular regions for ILM and RPE surfaces using the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 6 ; 6 3 ; 6 9 7 where Ω ILM 1.5 mm and Ω RPE 1.5 mm are the sets consisting of all faces within the 1.5-mm region of the ILM and RPE surfaces from their centers. The other parameters were defined in Eq. (9). Then, the total volume region is calculated using Eq. (16): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 7 ; 6 3 ; 6 0 1 Ω tv ¼ Ω 1.5 mm \ Ω cup ; (17) where Ω tv represents the total volume region on the ILM surface, as shown in Fig. 8(b), which is further employed for the total volume computation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 8 ; 6 3 ; 5 3 6 V tv ¼ where A i is the area of the face f i ∈ Ω tv and h i is the height with respect to corresponding vertex in the RPE surface similar to Eq. (12).

ONH annular region volume
The ONH annular region represents the ONH outer region, see Fig. 8(c). On ILM surface, this region is computed using the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 6 3 ; 3 9 6 Ω av ¼ Ω 1.5 mm \ Ω BMO ; (19) where Ω av consists of all the ILM surface faces, which belong to the annular region of the ONH. Similar to Eq. (18), we compute the volume of the annular region using the ILM and RPE surfaces correspondence: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 6 3 ; 3 2 1 V av ¼ The annular region volume helps to see the change in the outer region of the ONH volume in different cohorts.

Bending energy
The roughness on the ILM surface within the BMO region is an important parameter and commonly known as the bending energy on a manifold surface. The bending energy measures the fairness of a surface in terms of the curvature. In general, the outer region of the ILM surface is quite smooth and flat unlike the one inside the BMO, which has very complex topological structure. In this paper, we define the bending energy within the BMO region using the element-based normal voting tensor (ENVT). 44 The ENVT exploits the orientation information (face normals) to compute a shape analysis operator at each face f i ∈ Ω BMO and is defined as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 3 2 6 ; 6 1 3 where n j represents the normal of face f i and n T j is the transpose of n j . The term a j is the area of the face f i . To assure robustness against irregular sampling of the ILM surface, we weight Eq. (21) by the corresponding face area a j . The ENVT, M i , is a symmetric and positive semidefinite matrix, so, it can be decomposed into its spectral components: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 3 2 6 ; 4 9 9 where λ i k ¼ fλ i 0 ; λ i 1 ; λ i 2 g are the eigenvalues vector, and these eigenvalues are sorted in decreasing order ( . The corresponding eigenvector is denoted by e k . In general, the dominant eigenvalue λ i 0 has the corresponding eigenvector in the direction of the face normal and the remaining two eigenvectors will be aligned to the principle curvature direction on the ILM surface. On the planer region, only λ i 0 will be significant, on the edge region, λ i 0 and λ i 1 will be significant and at the corners, all of these eigenvalues are significant. Using the anisotropic properties of these eigenvalues, we define the bending energy inside the BMO region using the following equation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 3 ; 3 2 6 ; 3 2 8 where λ i 1 and λ i 2 are the two least dominant eigenvalues of the ENVT of the face f i . Figure 7(d) shows how each face of the BMO region is colored based on Eq. (26). The color is scaled from white (flat regions) to red (sharp features). We combine the ellipse fitted BMO points P 2D and the z-coordinates from P and represent these as P e ¼ fp e i ∈ R 3 ji ¼ 0; · · · ; n p − 1g, as it can be seen in Fig. 9 (green points). For each point p e i ∈ P e , we compute a point p MRW i on the ILM surface: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 ; 6 3 ; 5 4 2 p mrw The MRW points P MRW ¼ fp MRW i ∈ R 3 ji ¼ 0; · · · ; n p − 1g are lying on the ILM surface as shown in Fig. 9(a) (violet points). We create a quad surface using point sets P e and P mrw by introducing edges between the corresponding vertices in both point sets and connecting the neighbor points as shown in Fig. 9(b). The number of quad elements in the MRD surface is equal to the number of points in each point sets and is represented as Q ¼ fq i ji ¼ 0 · · · n p − 1g. MRW-MRA is computed as follows: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 6 ; 6 3 ; 4 0 9 A mrw ¼ where a q i represents the area of the quad q i .

BMO area
As shown in Fig. 5, BMOA represents the area under the fitted ellipse to the BMO points and is computed using the conic representation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 7 ; 6 3 ; 2 9 0 A BMO ¼ πr 1 r 2 ; where r 1 and r 2 are the major and minor axes of the fitted ellipse.

Experiments and Results
In order to evaluate our method, we performed repeated-measurement reliability tests, investigated ONH shape in healthy subjects, and tested if our method is able to detect differences in patients with diseases known to affect the ONH in the form of swelling and atrophy. Finally, we measured the implementation's performance.

Repeated-Measurement Reliability
In order to estimate the repeated-measurement reliability, we took three repeated scans of each eye from 10 healthy subjects. These subjects were measured each in a time frame of a week and then again in the following week. Table 1 shows the repeatability results. We see that our method scores highly at every parameter presented, with lowest intraclass correlation coefficient (ICC) of 0.905 for CONHT, and highest 0.998 for V cup . The ICC and confidence intervals were estimated using the variance components from a one-way ANOVA. We also evaluated our method with several other scan protocols of the same device (ONH cube with 73 B-scans, scanning angle of 15 deg ×15 deg and resolution 384 A-scans per B-scan, spatial resolution in x direction is ≈12.6 μm, in axial direction ≈3.9 μm and the distance between two B-scans ≈61 μm, ONH star scan with 24 B-scans, scanning angle of 15 deg ×15 deg and a resolution of 768 A-scans per B-scan, spatial resolution in x direction is ≈5.36 μm, in axial direction ≈3.9 μm) and the volumetric ONH-centered protocol acquired using Cirrus HD OCT (Carl Zeiss Meditec, Dublin, California), which covers 6 × 6 × 2 mm 3 region with 200 × 200 × 1024 voxels and obtained positive results.

Validation
We validated the BMO detection and checked the RPE segmentation. Five scans from the ones used in the repeatability testing were randomly selected. An experienced grader manually selected the BMO points using the RoI tool from ImageJ. 45 This resulted in a total amount of 488 B-scans with manually selected BMO points, which corresponded to the number detected automatically. Furthermore, we compared the mean signed and unsigned error in the x axis, as well as in the axial one (in z axis). If the automated method identified the BMO closer to the optic disc center, the sign of distance in the x-direction was positive. Similarly, if the automated BMO located below the manual BMO, the sign of distance in the z-direction was positive. Results are shown in Table 2. To our knowledge, there is no other study regarding BMO segmentation with the device used in this approach. As such, we are not able to relate our results to other published methods, since comparative studies revealed that measurements are not directly comparable between different OCT devices. 46 The ILM validation was done in another work of our research group. 40 For RPE computation, an experience grader visually controlled all scans used in the repeatability testing and assessed these as correctly segmented.

Clinical Evaluation
In this section, we present the results of our automatic pipeline approach for 248 OCT scans, from three groups, 71 HC eyes, 31 eyes of patients suffering from IIH, 25 which causes swelling of the ONH (papilledema). We also included 146 eyes of patients with autoimmune central nervous system disorders (MS and NMOSD) and a history of ON, an inflammatory optic neuropathy that damages the optic nerve leading to neuroaxonal degeneration. 47,48 All participants gave written informed consent according to the 1964 Declaration of Helsinki. The study was approved by the ethics committee of Charité -Universitätsmedizin Berlin.
In IIH patients, the ONH volume is increased and was shown to correlate with cerebrospinal fluid pressure. 8 The longitudinal analyses from Ref. 8 revealed that ONH volume measured by OCT decreased after the initial lumbar puncture and initiation of therapy with acetazolamide. Additionally, increased ONH volume was associated with lower visual acuity in IIH patients, which points out to the potential clinical relevance of the parameter.
ON is one of the most common initial clinical presentations of MS without any prior history of a demyelinating event. During the course of the disease, acute ON affects 50% and 70% of MS patients. 9 After initial swelling due to edema in the acute phase of ON, RNFL thickness decreased over the following six months. 47,49 ON is the first NMOSD-related clinical event in 55% of the patients, which causes severe structural damage to the optic nerve and retina with resulting functional impairment. Recurrent ONs in NMOSD give rise to severely thinned pRNFL and combined ganglion cell layer and inner plexiform layer. 11 Furthermore, within 5 years, ∼50% of NMO patients are either blind in one or both eyes. 49 RNFL reduction was consistently observed in MS patients, who had never had a clinical episode of ON, as well as in the clinically unaffected fellow eye of patients with a history of ON. 50,51 Similarly, a recent study 52 found that there is also retinal neuroaxonal damage without ON in NMOSD, but longitudinal studies investigating neuroaxonal damage without ON in NMOSD have not been performed extensively.
We hypothesized that patients with IIH would show significant changes in comparison to HCs indicating ONH swelling. Also, patients with a history of ON would show significant changes in comparison to HCs indication ONH atrophy. The results presented in Table 3 demonstrate that our approach successfully captures the differences. The only parameter showing no difference between groups is bending energy, E b . Although we expected that in the IIH affected eyes, the ONH shape inside the BMO to have a smoother convex shape, the data are still extremely heterogeneous. Thus, the bending energy reflects this extreme variability in the data. Especially inside the BMO, the topologies from one ONH to the other can be extremely different.
We plan to further investigate these findings (for both processes of swelling and atrophy of the ONH) with the parameters presented, to gain more insight into the pathology of the ONH.

Performance
Computation of the whole pipeline for a volume on an Intel (R) Core(TM) i7-4790, 3.60 GHz, 16 GB RAM and 64bit operating system takes on average 2 min (computation of the ILM, RPE lower boundary and BMO points takes about 75 s, and the morphological computation another 40 and 45 s).

Discussion and Conclusion
We have proposed an automatic and robust pipeline for the computation of several parameters that characterize ONH shape. In particular, we developed an approach that computes truly 3-D parameters derived from triangulated mesh surfaces of the ILM and RPE, and BMO points.
The proposed method identifies the RPE lower boundary and is able to provide a smooth 2.5-D surface by employing a twostage TPS fitting that preserves the retina natural curvature. Then, the estimated location of ONH region is detected and BMO points are identified. The issue of the presence of border tissue and the blood vessels was also dealt with by employing several filters that suppress vessel artifacts and enhance the Bruch's membrane termination points. Additionally, the computation time is significantly reduced by performing the BMO points' detection only in the ONH approximated region.  Furthermore, the ILM detection, using a Chan-Vese level set approach, provides a 3-D surface that is able to identify this membrane despite high heterogeneity and complicated topologies of the ONH. The issue of presence of noise and artifacts in the ILM surface, effects of acquisition process and marching cube algorithm, is corrected by employing a robust mesh denoising method. The resulting high fidelity surface preserves the features with minimal shrinkage. It is important to remember that we are interested in detecting even minimal shape changes; obtaining a surface that has no staircase artifacts, but preserves features, is crucial for parameter computation.
The same argument applies to the BMO points. We need a robust detection of this structure. Using our custom scan protocol, we have 145 B-scans in one volume, but only a part of these contain BMO points. Thus, depending on the position of each B-scan and on vessel shadow artifacts, large variation of BMO points in the XY plane can occur. To penalize the deviation of single BMO point coordinates from the elliptic shape of the BMO described in the literature, we fit an ellipse to the initially detected points. Hence, the natural trend of the BMO trajectories in the volume is kept while increasing the number of sample points.
Having these three components, ILM, RPE, and BMO points, our method computes the correspondence between surfaces ILM and RPE, such that each face of the ILM surface has a corresponding vertex in RPE. After this step, we have all the ingredients to derive the parameters.
Although the computed parameters, except the bending energy, were already introduced in several previous studies (mostly for the investigation of glaucoma or papilledema), we are the first to provide quantifiable measurements derived in a truly 3-D context. Furthermore, we were able to provide a first proof of applicability in a clinical setting, by showing the power of several parameters introduced in differentiating between inflammation in form of swelling and destruction in the form of reduction of layer thickness.
To the best of our knowledge, this is the first approach that is able to provide truly 3-D parameters. Most of the morphometric parameters derived for the characterization of the ONH shape 17,21,22,24 use either manual detection or manual correction of commercial software, a time-consuming task prone to subjectivity especially in volume scans consisting of a high number of B-scans. Similarly, methods 30,31,53 performed a manual detection of the BMO points and automatic segmentation of the ILM, BM, and choroid layers to compute surface correspondence between multiple eyes. Future work will include creation of population normative parameters and statistics using larger cohorts. The goal is to better detect and understand shape changes or differences correlated with neuroinflamatory diseases. Also, shapes from a large group of eyes will be made comparable by a common surface (a mean template); surface correspondence and anatomical changes at different time points will also be detected. Since we are interested in small, localized changes that could be clinically important, for example, as indicators of disease or treatment success, we will investigate the derivation of additional morphometric parameters.
There are several limitations. First, the OCT data were not acquired relative to the foveal-BMO axis, because some of our early data were taken before this concept had been established. OCT devices with wide field optics capable of imaging ONH and fovea in the same scan session may be used to address Table 3 Analysis of all the 3-D parameters defined for the HC and patient group. The last column shows the GEE analysis between the two groups.   (10) this issue in the future and establish the foveal-BMO axis as reference axis to eliminate rotational artifacts, which might increase repeated measurement reliability also in the context of longitudinal measurements. Furthermore, the algorithm was implemented and tested only on images of one device. Future work needs to adapt the approach to imaging data from other devices. Last, the sample size of our HC cohort is too small to allow meaningful investigation of ONH shape parameters correlations, e.g., with age and sex. These analyses should be repeated in a future study with higher sample sizes. We plan to perform an evaluation of our method on a larger dataset with different scan protocols and different OCT devices and to create a normative database for the developed parameters. In summary, our results showed that the proposed method successfully and robustly identifies ILM, RPE, and BMO points and it enables computing the morphometric parameters, that capture shape characteristics of the ONH.

Disclosures
A patent application has been submitted for the method described in this paper. Friedemann Paul serves on the scientific advisory board for the Novartis OCTIMS study; received speaker honoraria and travel funding from Bayer, Novartis, Biogen Idec, Teva, Sanofi-Aventis/Genzyme, Merck Serono, Alexion, Chugai, MedImmune, and Shire; is an academic editor for PLoS ONE; is an associate editor for Neurology Neuroimmunology and Neuroinflammation; consulted for Sanofi Genzyme, Biogen Idec, MedImmune, Shire, and Alexion; and received research support from Bayer, Novartis, Biogen Idec, Teva, Sanofi-Aventis/Genzyme, Alexion, and Merck Serono. Alexander U. Brandt served on the scientific advisory board for the Biogen Vision study; received travel funding and/or speaker honoraria from Novartis Pharma, Biogen, Bayer, and Teva; has consulted for Biogen, Nexus, Teva, and Motognosis, Nocturne. Sunil Kumar Yadav, Ella Maria Kadas, Seyedamirhusein Motamedi, Konrad Polthier, Frank Haußer, Kay Gawlik have no conflicts of interest to disclose.