1.1.

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.¹⁴ 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.¹⁷ and Muth and Hirnei¹⁸ 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.¹⁹ proposed BMO-MRW as a marker for early glaucoma detection also using manual segmented BMO points and BMO-MRW measurements, while Pollet-Villard²⁰ 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.²¹ analyzed BMO-MRW, the orientation of the long axis of BMO in a normal population, while Enders et al.¹⁷ described BMO-MRW in micro- and macrodiscs.²² 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²⁵ 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.²⁸ 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.²⁹ 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.³⁰ 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.³¹

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

1.2.

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:

2.1.

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\mathrm{deg}\times 15\mathrm{deg}$ and a resolution of 384 A-scans per B-scan. The spatial resolution in $x$ direction is $\approx 12.6\mu \mathrm{m}$, in axial direction $\approx 3.9\mu \mathrm{m}$, and the distance between two B-scans is $\approx 33.5\mu \mathrm{m}$.

2.2.

RPE Lower Boundary Surface Computation

2.2.1.

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 ($\sigma =5$ pixels isotropic with $\text{kernel size}=(10\mu \mathrm{m}\times 14\mu \mathrm{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³² 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.

Fig. 2

The first step in our pipeline, specifically the RPE lower boundary surface and BMO points detection exemplified using a B-scans. (a) One B-scan of the original volume (cyan arrows indicate blood vessels and the shadow artifacts these produce; the region delimited by the yellow line is part of the ONH disc). (b) After applying smoothing and intensity normalization. (c) Approximated ISOS junction points (blue points) after removing outliers. As a convention in our computation, the outliers detected at this step are set to have in the axial ($z$) direction the coordinate equal to 1. (d) Approximated upper boundary RPE points (blue points). (e) Smoothed 2-D RPE lower boundary (blue line). (f) BMO points (red dots); even in the presence of blood vessel the BMO points can still be detected. 384 pixels represent $\approx 4402.80\mu \mathrm{m}$ and 226 pixels $\approx 881.40\mu \mathrm{m}$.

2.2.4.

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, $\alpha$, 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 $\alpha =[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\times 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 $\alpha =[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 ${\mathrm{RPE}}_{\text{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 ${\mathrm{RPE}}_{\text{lower}}$ without including $A_{\mathrm{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).

Fig. 3

(a) RPE lower boundary surface with motion artifacts. (b) Resulting TPS fitted surface. The proposed two-stage TPS is capable of reducing the motion and segmentation artifacts seen in the original RPE lower boundary and creating a smoothed surface while still keeping the shape of the retina.

2.3.

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.¹⁴ 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,¹⁴ a structure similar to the BM. We propose a BMO points’ segmentation approach in the 3-D volume directly without the use of a 2-D projection image in the $xy$ plane, as presented in several previous works.^36,37

2.4.

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.⁴⁰ 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,⁴¹ where the input level sets are computed using the method proposed by Gawlik et al.⁴⁰ 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 ${\mathcal{M}}_{\mathrm{ILM}}$ with proper neighborhood and properly oriented face normals.

Fig. 4

(a) Noisy ILM surface. (b) Smooth ILM surface. (c) Smooth ILM and RPE surfaces with BMO. The ILM surface computed using the marching cube algorithm and level sets produced by the method.⁴⁰

2.5.

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.⁴² 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.

2.6.

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 $\mathcal{P}$ are projected onto the corresponding $XY$-plane and denoted as a points set: ${\mathcal{P}}_{2D}=\{{\tilde{\mathbf{p}}}_i\in {\mathbb{R}}^2|i=0,\cdots ,n_p-1\}$. Then, we apply the method⁴³ to compute a fitted ellipse to the BMO points in ${\mathbb{R}}^2$. Another key parameter in the ONH shape analysis is the center of the BMO points. This is computed as the barycentric of the all ${\mathcal{P}}_{2D}$ points:

Fig. 5

(a) Noisy and irregular BMO points, which are computed automatically. (b) Smooth and uniform density BMO points are computed using an ellipse fit.

2.7.

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\in {\mathbf{F}}_{\mathrm{ILM}}|i=0,\cdots ,m_{\mathrm{ILM}}-1$) of the ILM surface. In general, the RPE surface, represented here as a function graph: ${\mathcal{M}}_{\mathrm{RPE}}:{\mathbb{R}}^2\to \mathbb{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:

Fig. 6

Correspondence between ILM and RPE surfaces; the face $f_i\in {\mathbf{F}}_{\mathrm{ILM}}$ (green triangle) has a corresponding vertex ${\mathbf{v}}_i\in {\mathbf{V}}_{\mathrm{RPE}}$ in the RPE surface.

2.8.

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 ${\mathbb{R}}^2$ along with their center, as mentioned in Eq. (4). First, we compute the center of the ILM and the RPE surfaces corresponding to the center ${\tilde{\mathbf{p}}}_c$ of the BMO points. Let us consider that ${\tilde{\mathbf{V}}}_{\mathrm{ILM}}$ and ${\tilde{\mathbf{V}}}_{\mathrm{RPE}}$ represent projected sets of vertices onto $XY$-plane. Then, the indices of the centers of both surfaces can be computed as follows:

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 consider ${\tilde{\mathbf{p}}}_i^e$ represents the fitted ellipse point to ${\tilde{\mathbf{p}}}_i$ as mentioned in Sec. 2.6 and ${\tilde{\mathbf{p}}}_i^c$ denotes the corresponding affine transformed point on a circle of radius $r$. This transformation reduces the complexity of the BMO region computation and 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:

2.9.

3-D Shape Parameters

2.9.1.

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\in {\mathbf{F}}_{\mathrm{ILM}}$ and compute its centroid ${\mathbf{c}}_i$. As we discussed in Sec. 2.7, each face $f_i$ of the ILM has the corresponding vertex ${\mathbf{v}}_i\in \mathcal{C}$ in RPE. If $c_i^z-v_i^z\ge 0$ for all faces in ILM, then, there is no cup available in the ONH region. Otherwise, there is a cup. The terms $c_i^z$ and $v_i^z$ show the corresponding $z$-coordinates (height). Similarly, we can compute the cup region:

Eq. (10)

$${\mathrm{\Omega }}_{\mathrm{cup}}=\{f_i\in {\mathbf{F}}_{\mathrm{ILM}}|(c_i^z-v_i^z)\le 0\},$$

where ${\mathrm{\Omega }}_{\mathrm{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:

Eq. (11)

$$V_{\mathrm{cup}}=\sum _{(f_i\in {\mathrm{\Omega }}_{\mathrm{cup}})}{\mathbf{A}}_i{h}_i,$$

where ${\mathbf{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:

Eq. (12)

$${\mathbf{A}}_i=\frac{1}{2}({\overrightarrow{\mathbf{e}}}_0\times {\overrightarrow{\mathbf{e}}}_1),h_i=(c_i^z-v_i^z),$$

where ${\overrightarrow{\mathbf{e}}}_0$ and ${\overrightarrow{\mathbf{e}}}_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.

Fig. 7

A visual representation of the different morphological parameters. (a) BMO region, (b) cup volume $V_{\mathrm{cup}}$, (c) central ONH thickness (CONHT), and (d) bending energy $E_b$.

2.9.4.

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 ${\mathbf{v}}_{\text{ilmc}}$ and ${\mathbf{v}}_{\text{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:

Eq. (16)

$$\begin{gathered} {\mathrm{\Omega }}_{1.5\mathrm{mm}}^{\mathrm{ILM}}=\{f_j\in {\mathbf{F}}_{\mathrm{ILM}}||{\tilde{\mathbf{c}}}_j^{\mathrm{ILM}}-{\tilde{\mathbf{v}}}_{ilmc}|\le 1.5\mathrm{mm}\}, \\ {\mathrm{\Omega }}_{1.5\mathrm{mm}}^{\mathrm{RPE}}=\{f_j\in {\mathbf{F}}_{\mathrm{RPE}}||{\tilde{\mathbf{c}}}_j^{\mathrm{RPE}}-{\tilde{\mathbf{v}}}_{rpec}|\le 1.5\mathrm{mm}\}, \end{gathered}$$

where ${\mathrm{\Omega }}_{1.5\mathrm{mm}}^{\mathrm{ILM}}$ and ${\mathrm{\Omega }}_{1.5\mathrm{mm}}^{\mathrm{RPE}}$ 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):

Eq. (17)

$${\mathrm{\Omega }}_{tv}={\mathrm{\Omega }}_{1.5\mathrm{mm}}\backslash {\mathrm{\Omega }}_{\mathrm{cup}},$$

where ${\mathrm{\Omega }}_{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:

Eq. (18)

$$V_{tv}=\sum _{(f_i\in {\mathrm{\Omega }}_{tv})}{\mathbf{A}}_i{h}_i,$$

where ${\mathbf{A}}_i$ is the area of the face $f_i\in {\mathrm{\Omega }}_{tv}$ and $h_i$ is the height with respect to corresponding vertex in the RPE surface similar to Eq. (12).

Fig. 8

A visual representation of total and annular volume regions. (a) 1.5 mm regions of ILM (${\mathrm{\Omega }}_{1.5\mathrm{mm}}$) and RPE (${\mathrm{\Omega }}_{1.5\mathrm{mm}}^{\mathrm{RPE}}$) surfaces from the BMO center. (b) Total volume ($V_{tv}$) between ILM and RPE surfaces within 1.5 mm radius region. (c) Annular volume ($V_{av}$) is defined as the volume between ILM after taking out the BMO region.

2.9.8.

BMO-MRW surface area

BMO-MRW surface area, BMO-MRA, is computed by taking the whole region defined by all BMO-MRW. We combine the ellipse fitted BMO points ${\mathcal{P}}_{2D}$ and the $z$-coordinates from $\mathcal{P}$ and represent these as ${\mathcal{P}}_e=\{{\mathbf{p}}_i^e\in {\mathbb{R}}^3\hspace{0.17em}|i=0,\cdots ,n_p-1\}$, as it can be seen in Fig. 9 (green points). For each point ${\mathbf{p}}_i^e\in {\mathcal{P}}_e$, we compute a point ${\mathbf{p}}_i^{\mathrm{MRW}}$ on the ILM surface:

Eq. (25)

$${\mathbf{p}}_i^{mrw}={\mathbf{v}}_j\in {\mathbf{V}}_{ilm}|\min |{\mathbf{v}}_j-{\mathbf{p}}_i^e|.$$

Fig. 9

BMO-MRW surface with and without ILM. (a) BMO-MRW surface with ILM surface and (b) BMO-MRW surface.

The MRW points ${\mathcal{P}}_{\mathrm{MRW}}=\{{\mathbf{p}}_i^{\mathrm{MRW}}\in {\mathbb{R}}^3\hspace{0.17em}|i=0,\cdots ,n_p-1\}$ are lying on the ILM surface as shown in Fig. 9(a) (violet points). We create a quad surface using point sets ${\mathcal{P}}_e$ and ${\mathcal{P}}_{\text{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 $\mathcal{Q}=\{q_i|i=0\cdots n_p-1\}$. MRW-MRA is computed as follows:

Eq. (26)

$$A_{\text{mrw}}=\sum _{q_i\in \mathcal{Q}}{a}_{q_i},$$

where $a_{q_i}$ represents the area of the quad $q_i$.

3.1.

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_{\mathrm{cup}}$. The ICC and confidence intervals were estimated using the variance components from a one-way ANOVA.

Table 1

Repeatability test for the 3-D parameters.

We also evaluated our method with several other scan protocols of the same device (ONH cube with 73 B-scans, scanning angle of $15\mathrm{deg}\times 15\mathrm{deg}$ and resolution 384 A-scans per B-scan, spatial resolution in $x$ direction is $\approx 12.6\mu \mathrm{m}$, in axial direction $\approx 3.9\mu \mathrm{m}$ and the distance between two B-scans $\approx 61\mu \mathrm{m}$, ONH star scan with 24 B-scans, scanning angle of $15\mathrm{deg}\times 15\mathrm{deg}$ and a resolution of 768 A-scans per B-scan, spatial resolution in $x$ direction is $\approx 5.36\mu \mathrm{m}$, in axial direction $\approx 3.9\mu \mathrm{m}$) and the volumetric ONH-centered protocol acquired using Cirrus HD OCT (Carl Zeiss Meditec, Dublin, California), which covers $6\times 6\times 2{\mathrm{mm}}^3$ region with $200\times 200\times 1024$ voxels and obtained positive results.

3.2.

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.⁴⁵ 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.⁴⁶ The ILM validation was done in another work of our research group.⁴⁰ For RPE computation, an experience grader visually controlled all scans used in the repeatability testing and assessed these as correctly segmented.

Table 2

Mean unsigned and signed error in pixel and μm, for the x axis, and z axis between automatic (proposed) and manual segmentation.

3.3.

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,²⁵ 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.⁸ 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.⁹ 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.¹¹ Furthermore, within 5 years, $\sim 50\%$ of NMO patients are either blind in one or both eyes.⁴⁹ 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⁵² 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_{\mathrm{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.

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.

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.

3.4.

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