Layer boundary evolution method for macular OCT layer segmentation

Abstract

Optical coherence tomography (OCT) is used to produce high resolution depth images of the retina and is now the standard of care for in-vivo ophthalmological assessment. It is also increasingly being used for evaluation of neurological disorders such as multiple sclerosis (MS). Automatic segmentation methods identify the retinal layers of the macular cube providing consistent results without intra- and inter-rater variation and is faster than manual segmentation. In this paper, we propose a fast multi-layer macular OCT segmentation method based on a fast level set method. Our framework uses contours in an optimized approach specifically for OCT layer segmentation over the whole macular cube. Our algorithm takes boundary probability maps from a trained random forest and iteratively refines the prediction to subvoxel precision. Evaluation on both healthy and multiple sclerosis subjects shows that our method is statistically better than a state-of-the-art graph-based method.

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.

Fig. 2

A flowchart of our proposed method.

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.

Fig. 4

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

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.

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

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.

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.

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.

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.