1. Introduction

Optical coherence tomography (OCT) is a three-dimensional (3D) imaging technique that emits a beam of light into the tissues to be examined and detects the reflected or back-scattered light from the tissues. This reflected light will interfere with light that originated from the same source and the reflectivity profile along the light beam can be derived from the interference signal and used to generate an A-scan. Sequences of A-scans are used to form two-dimensional (B-scans) and three-dimensional volumes. The reflectivity profile of all A-scans at a fixed depth forms a C-scan. OCT acquires images with high resolution ($<5\mu$m) [1], from which the individual retina layers can be demarcated, qualitatively assessed, and quantified [2].

Ophthalmological assessment of retinal disease and neurological disorders often includes a segmentation of retinal cellular layers from OCT images [3–5]. Eight retinal subdivisions are commonly segmented: retinal nerve fiber (RNFL), ganglion cell and inner plexiform complex (GCIP), inner nuclear layer (INL), outer plexiform layer (OPL), outer nuclear layer (ONL), inner segment (IS), outer segment (OS), and the retinal pigment epithelium (RPE) complex. The segmentation of these eight layers are accomplished by identifying nine boundaries: vitreous-RNFL or internal limiting membrane (ILM), RNFL-GCIP, GCIP-INL, INL-OPL, OPL-ONL, ONL-IS or external limiting membrane (ELM), IS-OS, OS-RPE, and RPE-choroid or Bruch’s membrane (BM). As a result, an OCT image is partitioned into ten regions (see Fig. 1). We note that there is some debate as to the appropriateness of the term IS-OS, with some suggestion that ellipsoid zone (EZ) may be more correct [6]. Since this is an unresolved matter, we have continued to use the term IS-OS.

Quantification of layer thickness is often used in the study of various ophthalmological and neurological diseases. For example, subjects with multiple sclerosis (MS) of all subtypes exhibit significant thinning of the inner retinal layers [3]. OCT is also used to study myopia [7], diabetes [5, 8–10], Alzheimer’s disease (AD) [11], and glaucoma [4, 12, 13]. Since manual segmentation of each layer is time consuming and prone to inter-rater variability, fast, automatic algorithms for retinal layer segmentation are preferred.

 

Fig. 1 A B-scan from a healthy subject with manual delineation of nine boundaries. The boundaries from top to bottom are ILM, RNFL-GCIP, GCIP-INL, INL-ONL, OPL-ONL, ELM, IS-OS, OS-RPE, and BM.

Download Full Size | PDF

Various approaches have been developed for automated layer segmentation of OCT images [14–35], including methods based on shortest paths [17], active contours [18–20], statistical models [21], and level sets [22–25]. The graph method introduced by Li et al. [36], was originally developed for blood vessel segmentation before being repurposed by Garvin et al. [26] for OCT data. The underlying graph framework, with various improvements and extensions [14–16, 28–31], has become the predominate method in use today for OCT segmentation. Recently, a great amount of effort has been focused on exploring the use of deep networks, which have proven to be very effective in other medical imaging tasks, for use in OCT image segmentation. Fang et al. [33] developed a combined approach using both a deep neural network and a graph cut search, while He et al. [34] proposed a cascaded fully convolutional neural network to guarantee topological ordering of the layers.

 

Fig. 2 A flowchart of our proposed method.

Download Full Size | PDF

This paper introduces a novel framework for retinal layer segmentation that builds on the multi-layer fast level set method (MLS) [37]. Our method takes probability maps for each boundary from a trained random forest as input. An external force field and an initial boundary location are calculated from the input probability maps. The initialization is further refined through iterative updates based on the external forces and internal forces. Five major extensions to MLS are employed: 1) a subvoxel location is calculated during surface evolution, which allows for better estimation of the smoothness force; 2) a simple Gaussian filtering process is used to calculate the smoothness force without the need to calculate curvature, which significantly reduces computational time; 3) the initial surface is generated using a shortest path method; 4) surface evolution is performed in a BM flattened space and the benefit of using a flat space is incorporated into an extra term in the internal force; and 5) vector field convolution [38] is used instead of gradient vector flow [39] for better performance. Our proposed layer boundary evolution (LBE) algorithm can be divided into two parts. The first part, in which nine boundary probability maps are calculated from a retinal OCT image, is largely adopted from [28]. Three steps are involved: preprocessing an OCT image, feature extraction, and boundary classification using a trained random forest. The second part uses the boundary probability maps from the trained random forest for generating two inputs to the layer boundary evolution: a force field and an initial prediction for each boundary. The final result is the output of layer boundary evolution, which iteratively refines the input prediction to a subvoxel accuracy (see Fig. 2). Since the first part is well documented in [28] and its source code is publicly available (http://www.nitrc.org/projects/aura_tools), we briefly introduce each step in the first part and then focus most of our discussion on the second part.

2. Method

2.1. Preprocessing

The preprocessing steps serve to transform an OCT volume into a common space so that the intensity and spatial coordinates of an input volume is consistent with the training data. An example of a preprocessed B-scan is shown in Fig. 3.

2.1.1. Retina boundary flattening

Flattening is a common preprocessing step performed by many retina segmentation algorithms [17, 26]. We first estimate the top and bottom boundaries of the retina, specifically the ILM and the BM. A B-scan is filtered by a Gaussian smoothing kernel with σ = 3 pixels and the derivative along each A-scan within the B-scan is computed. The locations of the ILM and BM in a particular A-scan is found as the locations that have the largest positive and negative derivatives, respectively. Error points are found by comparing the estimated boundaries with its independently median filtered version. Points that deviate more than 15 voxels from the median filtered surface are defined as outliers and are replaced by interpolated points. A BM flattened volume is obtained by vertically shifting each A-scan within its B-scan based on the estimated BM surface. Features related to spatial coordinates of pixels contain more information in a retina boundary flattened volume. Our classifier described in Section 2.2 is trained using features calculated from this flattened volume and its corresponding manual delineation.

 

Fig. 3 The outputs of our preprocessing steps. A B-scan from a healthy subject is shown (a). The ILM and BM are detected (b) and used to produce a BM flattened B-scan (c). The image in (d) shows a normalized flattened B-scan.

Download Full Size | PDF

2.1.2. Intensity normalization

The automatic intensity rescaling and automatic real-time averaging performed by the OCT scanner may potentially cause a large intensity variation for a particular tissue type within one subject. These differences in dynamic range are not desirable especially during feature classification. An intensity consistent B-scan is achieved by computing a robust maximum I_m and linearly rescaling intensities from the range $\left[0,I_m\right]$ to the range $\left[0,1\right]$, with values greater than I_m being clamped at 1. The robust maximum I_m is found by first median filtering each individual A-scan within the same B-scan using a kernel size of 15 pixels. Then, I_m is set to the value that is 5% larger than the maximum intensity of the entire median-filtered image.

2.2. Random forest classification

The boundary locations are estimated using a trained random forest classifier [28]. In total 27 features are calculated for every voxel in the preprocessed OCT volume. Those features are passed to decision trees in the forest, which provides information about whether the voxel belongs to a specific boundary or not. By combining the predictions from all decision trees, the probability of each voxel belonging to each boundary is computed. In short, the trained random forest takes features and outputs nine probability maps, each associated with one boundary (see Fig. 4). A probability map has the same size as the input OCT volume. To train the random forest, manual delineations are also preprocessed to have the same volume size and flattened BM.

 

Fig. 4 Shown are the random forest predicted boundary locations in a B-scan.

Download Full Size | PDF

2.3. Surface representation

The fast level set method introduced by Shi et al. [40] uses a list instead of a level set function to represent a boundary or the zero level set. Their idea adapts naturally to boundaries in retinal OCT images. Specifically, the anatomical structure of human retina determines that retinal layers cannot fold on themselves. In other words, an A-scan can only intersect a specific layer boundary once in an OCT image. Therefore, for a given OCT image of size $l\times m\times n$, where l is the number of C-scans, m is the number of A-scans within a B-scan, and n is the number of B-scans within the 3D volume, a 2D array of size $m\times n$ is sufficient to represent a specific boundary. The j^th row, k^th column of the 2D array corresponds to the j^th A-scan within the k^th B-scan, and its value indicates the intersection location of that boundary in that A-scan. In other words, this 2D array represents the boundary height as shown in Fig. 5. We note that a layer boundary in a B-scan of size $l\times m$ is represented by a 1D array of size $1\times m$.

 

Fig. 5 Representation of the ILM surface. The red line segment shows the height of the ILM (320 voxels in this case) at the ${350}^{\text{th}}$ A-scan of the ${40}^{\text{th}}$ B-scan. Therefore in the 2D array representation of the ILM surface, the value at location $\left(350,40\right)$ is 320. For a $496\times 1024\times 49$ OCT volume, since there are 1024 × 49 A-scans, the 2D array for the ILM surface is of size 1024 × 49. In total nine 1024 × 49 2D arrays are needed to fully represent all nine boundaries.

Download Full Size | PDF

2.4. Force field

The external force F_ext used to drive the surface to the layer boundaries is computed from the probability map described in Section 2.2. A probability map can be thought of as an edge map of the corresponding boundary. We generate the external force field using vector field convolution (VFC) [38], which calculates the convolution of the probability map and a vector field kernel. Typically, all vectors in the vector field kernel point to the kernel origin; thus our external force always points to locations with the highest boundary probability. It can be used to guide our model to the boundary locations learned by the random forest. An example of the vector field kernel with its generated force field in a zoomed B-scan is shown in Fig. 6. During layer boundary evolution, we evolve surfaces vertically in a B-scan (see Section 2.6), which means only forces along an A-scan take effect. Which is why we only show vertical forces. The vector field kernel, however, is 2-dimensional so that information from adjacent A-scans is incorporated. Compared with other methods, VFC produces a force field that is robust to noise. This is important for processing fuzzy edge maps, like the probability maps generated by a random forest. A comparison between VFC and gradient vector flow [39], which our previous approach used, found a three fold improvement in speed [38].

 

Fig. 6 The vector field kernel used is shown in (a). A probability map of ILM (b top) is convolved with this vector field kernel to produce the external force field. A magnified view of the force field (red arrows) is overlaid on the probability map and shown in (b bottom).

Download Full Size | PDF

Internal forces are used to enforce smoothness of the output surface. The smoothness constraints for layer boundaries in retinal OCT are different from other boundaries. Specifically, boundaries have very different curvature near the fovea versus far from the fovea. Because of this, a universal smoothness constraint over the entire boundary may produce overly-smooth surfaces near the fovea and bumpy surface elsewhere. Using a flat space [25, 41, 42] can solve this problem by encoding the boundary shape prior information within the flat space. A flat surface in the flat space corresponds to a curvy surface near the fovea and smooth surface away from the fovea. Therefore enforcing a universal smoothness constraint in the flat space allows different smoothness characteristics in different regions. However, a smooth surface in flat space does not necessarily mean a smooth surface in native space, because the conversion from flat space to native space involves interpolation [42]. Instead of converting to flat space, we propose a new internal force $F_{\text{int}}=F_{\mathrm{s}}+F_{\mathrm{p}}$ that accounts for the smoothness of the surface and uses a shape prior for the retinal boundaries.

Internal forces for smoothness regularization, F_s, are commonly computed using the curvature on the surface, which is computationally expensive to evaluate. Alternatively, we introduce a Gaussian filtering process that efficiently calculates F_s. Use of Gaussian smoothing to replace the curvature was first proposed in [40]. It is motivated by the observation that when the level set function is chosen as the signed distance function, the curvature takes a simpler form and equals the Laplacian of ϕ, From the theory of scale-space ([43]), we know that evolution of a function according to its Laplacian is equivalent to Gaussian filtering the function. We first filter the predicted surface B, where B is the 2D array of numbers representing the retinal surface as depicted in Fig. 5, with a Gaussian kernel of standard deviation s to produce a smoothed surface B_s. F_s is then computed as,

The second term F_p is computed using the difference between the smoothed initialization, B_p, and the predicted surface B as,

The external force F_ext computed from VFC is normalized between $\left[0,1\right]$. During the evolution described in Section 2.6, we compute the internal force F_int for a boundary based on the current predicted surface B. The total force F for evolution is computed as:

2.5. Initialization

With an external force field free of noise, generating an initialization is a trivial problem as the forces point perfectly towards the true boundary. For example, we could simply initial nine flat planes as the predicted surfaces. However, as a real force field is usually noisy, pixels away from the boundary may not have forces that point in the correct direction. If we initialize the surface far from the true boundary, it is likely that it will never evolve to the true boundary, because false boundary responses from noise in the neighborhood may overwhelm the force for a true boundary. Therefore, we want to initialize the surface close enough to the true boundary to guarantee that our initialization is within the capture range of the force fields. We use Dijkstra’s shortest path (DSP) algorithm [44] to find the shortest path using the B-scan probabilities.

To use DSP, we construct a graph within each B-scan for each boundary using its probability $\mathcal{P}\left(i,j\right)$ derived from the random forest. The vertex at the i^th location in the j^th A-scan only connects to vertices in the ${(j+1)}^{\text{th}}$ A-scan that are between $i-1$ and $i+1$. This ensures that the shortest path we find is smooth and, more importantly, it goes through each A-scan only once. The edge cost $\mathcal{C}\left(i,j\right)$ to the vertex at location i in the j^th A-scan is defined as

The total distance is the sum of all edge costs from the previous vertices. We initialize all vertices in the first A-scan as visited with the total distance equal to the edge cost described above. We then iteratively assess a vertex with the minimal total distance within all visited vertices. The algorithm stops when we mark a vertex in the last A-scan as visited. From that vertex, we trace back to the first A-scan. The coordinates of all vertices along the path are used as an initial surface location. For an input B-scan probability map, we find a smooth path that goes across all A-scans. Since the shortest path method provides no topology guarantee, its predicted surfaces may cross over each other, a problem most likely to occur near the fovea. To address this problem, after the nine boundary predictions are individually found in each B-scan, we fix ordering errors as we combine them. We first combine the boundaries of ILM and RNFL-GCL, if in some A-scans the boundary height of RNFL-GCL are close to or even larger then the boundary height of the ILM, we move the RNFL-GCL boundary in those A-scans down so that the RNFL layer has a predefined minimal thickness. This process is iterated to combine all nine boundaries, guaranteeing that the initial boundaries have the correct layer ordering along every A-scan. Previously, we used a regression model to generate initial boundary locations from the retina thickness [42]. We find that the DSP is fast and robust and produces an initialization that is closer to the true boundary.

2.6. Layer boundary evolution

Layer boundary evolution (LBE) is an iterative process that refines the initial surfaces to the boundary locations indicated by the force fields described in Section 2.4. Our evolution scheme evolves the boundaries in the A-scan direction by adjusting values in the 2D arrays. Adding to or subtracting from a 2D array will move the boundary towards the inner or outer retina, respectively. This approach provides two benefits: firstly, we can easily find grid points adjacent to the boundary without the need of a level set function; secondly, the layer ordering can be easily maintained during the evolution.

Given the layered nature of the tissues under consideration within the retina, there is a topological requirement for OCT images. Specifically, layers cannot fold over themselves and the ordering of the eight layers between the vitreous and choroid needs to be maintained. Since only one location is recorded in each A-scan for each boundary using our surface representation, the first requirement is automatically satisfied. We postprocess the initialization to have the correct layer ordering, as described in Section 2.5. Evolving the boundary only in the A-scan direction helps maintain the ordering of layers to be consistent with the underlying tissue ordering. Every time a change of boundary height (a value changes in the surface 2D array) that will result in its adjacent layer to have thinner thickness than the minimal thickness, we reset its boundary height so that layer has the exact minimal thickness defined. As is shown in Fig. 7, before the current iteration, Point C and Point D are $0.9$ pixels away. If Point C tries to move down by more than $0.8$ pixels, its location will be reset so that the distance between Point C and Point D is exactly $0.1$ pixels, where $0.1$ pixels is the predefined minimal layer thickness (for this example).

 

Fig. 7 Shown are two surfaces in front of a B-scan. We highlight two A-scans (the ${250}^{\text{th}}$ and the ${500}^{\text{th}}$ A-scans at the ${21}^{\text{th}}$ B-scan). The first A-scan intersects with the IPL-INL surface and INL-OPL surface at Points C and D; the second A-scan intersects with the IPL-INL surface and INL-OPL surface at Points A and B.

Download Full Size | PDF

During each iteration, a decision will be made for every value in a 2D array B: whether to increase the value (move up towards inner retina), decrease the value (move down towards outer retina), or keep the surface position unchanged. For a boundary point $B^n\left(j,k\right)$, this decision is made by looking at the combination of external and internal forces at its upper adjacent grid point ($\text{ceil}\left(B^n\left(j,k\right)\right)$) which is computed as

Function $\text{ceil}(\cdot )$ rounds the input value to the nearest larger integer while function $\text{floor}(\cdot )$ rounds the input to the nearest smaller integer. Superscript n indicates the iteration index, and $\left(j,k\right)$ is the A-scan index. While F_ext can be precomputed, we need to calculate F_int during each iteration. We use two arrows and a division line to represent the direction of the net forces $F_{\mathrm{u}}^n\left(j,k\right)$ and $F_{\mathrm{l}}^n\left(j,k\right)$; for example, the symbol $\frac{\Downarrow }{\Uparrow }$ indicates that the $F_{\mathrm{u}}^n\left(j,k\right)$ is pointing down while the $F_{\mathrm{l}}^n\left(j,k\right)$ is pointing up. The decision rule for refining a predicted surface point $B^n\left(j,k\right)$ in the n^th iteration is described as follows:

Case 1 has two net forces pointing at each other, which is the feature exhibited at a boundary. The zero level set location is calculated using the net force at the upper grid point $F_{\mathrm{u}}^n\left(j,k\right)$ and lower grid point $F_{\mathrm{l}}^n\left(j,k\right)$:

Even if this subvoxel location may change in future iterations because of changes in the internal force, it is essential to calculate it so that boundary points in adjacent A-scans can use this more accurate location for better estimating F_s in internal force in later iterations. For cases 2 and 3, we evolve the boundary point to follow the force field. Both cases indicating that the true boundary is not in between the current upper grid point and lower grid point interval. Situation 4 accounts for the influence of noise or a bad initial surface. When we have $F_{\mathrm{u}}^n\left(j,k\right)$ and $F_{\mathrm{l}}^n\left(j,k\right)$ pointing away from each other, the combined force field cannot correctly infer the direction of the true boundary. It can result from noise in the probability map or a poor initial surface. In such cases, our model trusts the smoothness force F_s. Because the upper grid point ($\text{ceil}\left(B\left(j,k\right)\right)$) is always one voxel higher then the lower grid point ($\text{floor}\left(B\left(j,k\right)\right)$), from Equation 1 it is obvious that the smoothness force at the upper grid point will always be smaller than the lower grid point. Thus, in subsequent iterations we will not have $F_{\mathrm{u}}^n\left(j,k\right)$ and $F_{\mathrm{l}}^n\left(j,k\right)$ pointing away from each other again.

3. Experiments and results

3.1. Data set

We test our LBE algorithm on macular OCT images with no observable intraretinal pathologies. We do this with data from the right eyes of 36 subjects obtained using a Spectralis OCT system (Heidelberg Engineering, Heidelberg, Germany). The research protocol was approved by the local Institutional Review Board, and written informed consent was obtained from all participants. Of the 36

subjects, 21 were diagnosed with MS while the remaining 15 were healthy controls. All scans were screened and found to be free of microcystic macular edema (a pathology sometimes found in MS subjects).

All scans were acquired using the Spectralis scanner’s automatic real-time function in order to improve image quality by averaging at least 12 images of the same location (and the ART setting was fixed at 12). The resulting scans had signal-to-noise ratios (SNRs) of between 20 dB and 38 dB, with a mean SNR of $30.6$ dB. Macular raster scans $\left({20}^{\circ }\times {20}^{\circ }\right)$ were acquired with 49 B-scans, each B-scan having 1024 A-scans with 496 voxels per A-scan. The B-scan resolutions varied slightly between subjects and averaged $5.8\mu$m laterally and $3.9\mu$m axially. The through-plane distance (slice separation) averaged $123.6\mu$m between images, resulting in an imaging area of approximately 6 × 6 mm². We note that our data is of a higher resolution and better SNR than many other recent publications. The nine layer boundaries were manually delineated on all B-scans for all subjects by a single rater using an internally developed protocol and software tool. The manual delineations were performed by clicking on approximately 20–50 points along each layer border followed by interpolation between the points using a cubic B-spline. Visual feedback was used to move points to ensure a correct boundary.

3.2. Parameter selection

The random forest classifier was trained using 4 healthy subjects and 6 MS subjects. The recommended parameters of the RF classifier described in [28] are used. Since our initialization is very close to the final results, the maximum number of iterations is set to 10. This allows an accurate estimation of the boundary at a small computational cost.

There are three parameters related to the internal force that need to be determined in our proposed method (LBE) for a single boundary:

Different sets of parameters were found for each of the nine boundaries using a validation set of 2 healthy subjects and 3 MS subjects. We used mean absolute boundary error as the metric when determining these parameters. The mean absolute boundary error is calculated as the absolute distance for each boundary against the manual delineated surface along each A-scan. These errors are first averaged across all A-scans in each subject and then averaged over all subjects for each boundary.

 

Fig. 8 Shown are the effect of parameter selection on absolute boundary error for four boundaries in the validation set. Each sphere corresponds to an experiment on a validation set and the parameters used are indicated by the spatial coordinates of the sphere. After all 1440 experiments, the resultant absolute error were thresholded and errors that were larger than 5% of the lowest error are shown as big red spheres. The remaining results were linearly mapped to $\left[0,1\right]$ with the lowest error indicated by 0 (small blue spheres) and largest errors indicated by 1 (large yellow spheres). The parameters chosen are indicated by a magenta triangle. The shaded green region indicates the combination of parameters that our method produce lower absolute boundary error than AURA v3.4 when evaluated on the test data set. The black contour indicates the intersection of this region and each slice.

Download Full Size | PDF

We ran 1440 experiments on the validation set for each boundary with different parameters and we chose our parameters by selecting the combination of c_s, c_p, and s that produced the lowest mean absolute boundary error. We limited the search region of c_s and c_p to be within $\left[0,2.75\right]$; otherwise, the internal force might have overwhelmed the external force yielding a result that is simply a smoothed version of the initial surface. Figure 8 shows the effect of parameter selection on the mean absolute error of four boundaries. Each point in the 3D volume in Fig. 8 represents the resultant absolute error produced by a certain choice of the three parameters. For each boundary, the resulting mean absolute errors larger than 5% of the lowest mean absolute error are shown as red circles. Meanwhile a mean absolute error within the range of 5% of the best result is linearly mapped to $\left[0,1\right]$ with the lowest mean absolute indicated by 0 (blue small spheres) and the worst mean absolute error indicated by 1 (yellow large spheres). The best results are marked with a magenta triangle. Figure 8 also indicates the region (green) of parameters where LBE produce a lower mean absolute error than a graph method (AURA v3.4) we compare with in Section 3.3 over the remaining 21 subjects. The black contour is the intersection of this region and each slice. The RNFL-GCIP boundary (top right subplot) has no green region indicating that no combination of the three parameters can produce better result than AURA v3.4, while for IS-OS boundary (bottom left subplot), any combination of parameters can produce better result than AURA v3.4. The figure for GCIP-INL and INL-OPL boundaries are very similar to ILM (top left subplot), the figure for OPL-ONL, ELM, and BM are similar to IS-OS (bottom left subplot), for space consideration the figures for these five boundaries are not shown.

After determining the optimal c_s, c_p, and s for the validation set, we used them on the test data set; the mean absolute boundary error for iterations between 1 and 50 are shown in Fig. 9. This demonstrates that 10 iterations is an appropriate number to allow all boundaries to converge.

 

Fig. 9 Shown is the mean absolute boundary error (in pixels) for all nine boundaries from iteration 1 to iteration 50. Each pixel has a physical dimension of $3.9\mu$m axially and $5.8\mu$m laterally.

Download Full Size | PDF

3.3. Results

We evaluated our method (LBE) on the remaining 21 subjects that were not used to train the random forest or for parameter selection. We compared with a previously reported graph cut method (AURA v2.11) [28] and its most recent version (AURA v3.4). AURA v2.11 has been independently shown to be one of the most accurate OCT retinal segmentation tools [45]. For the purpose of direct comparison, the recommended parameters [28] for the random forest are used for all three methods. An example of the manual delineations, as well as the result of the AURA v3.4 and LBE methods on the same subject are shown in Fig. 10.

 

Fig. 10 A B-scan retinal OCT (a) and a magnified region near the fovea (b) are shown, along with corresponding manual delineation (c) and (d). The B-scan is segmented using AURA v3.4 (e) and (f), and our method (LBE) (g) and (h). The LBE mask is generated by rounding the boundary location to $0.1$ voxel level.

Download Full Size | PDF

Shown in Table 1 is the mean absolute boundary error of the three methods against the manual delineation. The mean absolute boundary error is calculated as the absolute distance for each boundary against the manual delineated surface along each A-scan. The standard deviations were calculated assuming that every error value is separate (i.e., errors were not averaged for each subject before calculation). It is observed that the mean absolute boundary error is smaller for LBE than AURA v3.4 in 8 out of the 9 boundaries, the exception being the RNFL-GCIP boundary. A paired Wilcoxon signed rank test was used to test the significance of any improvement between AURA v3.4 and our method. 6 of the 9 boundaries reach significance (α level of $0.001$). Two boundaries (OS-RPE and BM) lack statistical significance because of the large variance. AURA v3.4 has lower mean absolute boundary error than LBE on one of the nine boundaries, however it does not reach significance.

Table 1. Mean and standard deviations of absolute boundary error (in pixels) are shown. Each pixel correspond a physical dimension of $3.9$ microns axially and $5.8$ laterally. A paired Wilcoxon signed rank test was used to test the significance of any improvement between AURA v3.4 and our method (LBE). The method with the lowest mean absolute boundary error are bolded. Boundaries that are significant (α level of $0.001$) are marked with a $\star$. We include the result directly from shortest path initialization (SP) to show that our method improves on the initialization. A paired Wilcoxon signed rank test (not shown in the table) between LBE and SP shows all boundaries except BM reach a p-value of 0.0001, where the p-value of BM is 0.0143 due to the large variance in the LBE output.

View Table | View all tables in this article

We also computed the Dice Coefficient [46] between each of the three methods and the manual delineation; the results are shown in Table 2. The Dice Coefficient measures how much the segmentation results produced by an algorithm agree with the manual delineation. The Dice Coefficient has a range of $\left[0,1\right]$, with higher values reflecting more agreement between the algorithm and the manual delineation. Since both manual delineation and segmentation results are boundary locations, we first converted the boundary locations to a layer mask by rounding the boundary location to the voxel level and then in each A-scan assigning the corresponding layer label from the current boundary voxel to the location for the boundary below. The resulting Dice coefficients for AURA v3.4 and LBE are very close, with five out of the eight layers showing statistically significant improvements (α level of $0.001$).

Table 2. Mean and standard deviations of the Dice Coefficient across the eight retinal layers. A paired Wilcoxon signed rank test was used to test the significance of any improvement between AURA v3.4 and our method (LBE). The method with the highest Dice Coefficient are bolded. Boundaries that are significant (α level of $0.001$) are marked with a $\star$.

View Table | View all tables in this article

The computational performance are compared between AURA v3.4 and LBE. Both algorithms were implemented in MATLAB 2018a and run on a four core 3.1GHz Linux computer. For our 21 test data set, the average computational time per subject for LBE is $93.4\left(\pm 7.4\right)$ seconds compared with $106\left(\pm 4.8\right)$ seconds for AURA v3.4. Typically, more than 80 seconds are spent on producing the probability maps, which are common to both algorithms. The peak memory usage for our algorithm is $9.7$ GB during the random forest classification, compared with $11.5$ GB for AURA v3.4 during its graph cut search.

4. Conclusion

In this paper, we proposed a layer boundary evolution method that loosely originates from the fast level set method to replace the commonly used graph based method. We build our pipeline upon the state-of-the-art graph based method [28]: our method takes boundary probability map from a trained random forest, the result produced by the proposed method shows excellent performance in terms of both mean absolute boundary error and layer Dice coefficient. In particular, five out of nine boundary estimates are statistically better than the state-of-the-art graph cut method when using the same random forest classifier. In addition, less computation cost is required, which makes it an ideal alternative for graph cut search in boundary refinement after feature classification. The mean absolute boundary error is less then $4\mu$m for all nine boundaries. Although we are unable to outperform deep learning based methods [47] in terms of speed (since the latter needs only a few seconds to run on a GPU), our method provides topological correct results and subvoxel resolution that most deep learning based methods cannot achieve. In terms of the run times of other methods, we note that OCTRIMA3D [48] has a run time of 28 seconds whereas AURA v2.11 on the same hardware had a run time of 152 seconds [48]. However, AURA v2.11 performed better in terms of mean absolute boundary error.

The data used to develop LBE is of higher quality than is typically acquired in a clinical OCT setting, however it is similar to the data used to developi the state-of-the-art [28]. Moreover, the random forest that generates the boundary classification probability maps is the same as [28] and has been previously shown to be very robust across scanners and data sources [45, 49]. Which provides us with some confidence about the general utility of out method across scanning platforms.

Future work will focus on replacing the current random forest classifier with a deep learning model to produce boundary probability maps using data driven features. The major challenge will be producing a diffused probability map from a deep network. Deep networks are often quite certain about the position of a boundary, with near binary probability maps as output; that is the true uncertainty of the deep network is not fully reflected in the networks output.