3D Automatic Segmentation Method for Retinal Optical Coherence Tomography Volume Data Using Boundary Surface Enhancement

With the introduction of spectral-domain optical coherence tomography (SDOCT), much larger image datasets are routinely acquired compared to what was possible using the previous generation of time-domain OCT. Thus, there is a critical need for the development of 3D segmentation methods for processing these data. We present here a novel 3D automatic segmentation method for retinal OCT volume data. Briefly, to segment a boundary surface, two OCT volume datasets are obtained by using a 3D smoothing filter and a 3D differential filter. Their linear combination is then calculated to generate new volume data with an enhanced boundary surface, where pixel intensity, boundary position information, and intensity changes on both sides of the boundary surface are used simultaneously. Next, preliminary discrete boundary points are detected from the A-Scans of the volume data. Finally, surface smoothness constraints and a dynamic threshold are applied to obtain a smoothed boundary surface by correcting a small number of error points. Our method can extract retinal layer boundary surfaces sequentially with a decreasing search region of volume data. We performed automatic segmentation on eight human OCT volume datasets acquired from a commercial Spectralis OCT system, where each volume of data consisted of 97 OCT images with a resolution of 496 512; experimental results show that this method can accurately segment seven layer boundary surfaces in normal as well as some abnormal eyes.

valleys in an A-scan intensity profile by using a mean filter for de-speckling, which segmented five layer boundaries [3]; Chiu et al. presented a segmentation method that used graph theory and dynamic programming to segment seven retinal layers [4];This method was later extended for segmentation of mouse retinal layers [5], anterior eye images [6] and age-related macular degeneration(AMD) images [7]; Using a similar method, Yang et al. [8,9] utilized a more complex approach to calculate the weights map of graph-based method, using dual-scale gradient information and shortest path search techniques to segment intra-retinal boundaries in OCT images. Yazdanpanah et al. used an active contour approach for the segmentation of rodent retinas [10]; A two-step kernel-based optimization was proposed by Mishra et al. [11]. However, the methods in [10,11] was never tested on OCT datadets of human retinas. Itebeddine et al. proposed a global segmentation algorithm based on using active contours and Markov random fields to segment eight retinal layers [12].
While the aforementioned methods can be used to segment OCT volume data slice-by-slice, most of them require long processing times, and they do not use the correlation between slices well. Recently, 3D OCT retinal image segmentation techniques have been developed. Zawadzki and Fuller et al. proposed segmentation methods using a support vector machine (SVM) and machine learning [13,14], which could segment one layer once by manual interaction. Similar to Zawadzki et al.' method, a support vector machine was used to classify pixels in the OCT image but in a fully automated way [15].
In [16], random forest classifier was built to segment eight retinal layers in macular cube images acquired by OCT. The random forest classifier learns the boundary pixels between layers, producing an accurate probability map for each boundary, which is then processed to finalize the boundaries. Kajic et al. proposed a method that used a large training dataset obtained from manual segmentations by human operators as input to develop a statistical model to segment seven retinal layers [17]. Garvin et al. proposed a graph search-based three-dimensional OCT retinal image segmentation algorithm [18,19], which could segment five retinal layers, which was later extended to incorporate hard/soft constraints [20]. Lee used multi-scale 3D graph search for segmenting the optic nerve head [21]. Compared to these complex 3D imaging segmentation approaches, Fabritius et al. presented a fast segmentation method for segmenting the internal limiting membrane (ILM) and the retinal pigment epithelium (RPE) that was based on variations in pixel intensity [22];While the method is relatively fast, only two boundaries are detected.
To develop a more practical 3D segmentation technique such as computation efficiency and robust to blood vessel shadow and noise, here we propose a new 3D segmentation method for segmenting retinal layer boundaries from OCT volume data using boundary surface enhancement, which is relative simple, efficient, robust to shadows and noise, and can segment seven boundary surfaces. To segment a boundary surface, two OCT volume datasets are obtained by using a 3D smoothing filter and a 3D differential filter. Their weighted sum is then calculated to generate new volume data with an enhanced boundary surface, where the pixel intensity, boundary position information, and intensity changes on both sides of the boundary surface are used simultaneously. Then, preliminary discrete boundary points are detected from the A-Scans of the volume data. Finally, surface smoothness constraints and a dynamic threshold are applied to obtain a smoothed boundary surface by correcting a small number of error points. Our methods can extract retinal boundary surfaces sequentially within a decreasing region of volume data. The key idea is to use pixel position information, gradient information and intensity information simultaneously to enhance the boundary surface to be detected so that preliminary discrete boundary points can be detected more correctly and error points can be eliminated more easily. This paper is organized as follows: Section 2 gives a description of our generalized layer segmentation algorithm, which is fundamental for segmenting all of the layer boundary surfaces; Section 3 demonstrates how to segment seven retinal layer boundary surfaces in detail; experimental results and analysis are given in Section 4; and conclusions are made in Section 5.

2． A Generalized Layer Segmentation Algorithm
An image acquired from a commercial Spectralis OCT device (Heidelberg Engineering, Heidelberg, Germany ) is shown in Fig. 1, where the left panel shows the scanning position in the retinal tissue and the right panel shows the corresponding SDOCT image. Fig. 2 illustrates a volume dataset made up of a sequence of SDOCT B-scans in the y direction. Every B-scan is composed of A-scans in the x direction. Each A-scan has a depth coordinate z, which increases going from top to bottom in the image.   This section proposes a method for segmenting the layer structures of the retina. The basic idea of this method is to use the characteristics of the boundary of interest to design a 3D operator and apply it to the original volume data so that the pixel value on the desired boundary in the new volume data is likely to be the maximum value in its A-scan.
This approach makes the new volume data a better indicator of the desired boundaries.
The algorithm consists of three steps: denoising, extracting boundary points and correcting error points. Obviously, the key problem is to determine how to identify the correct discrete boundary points. Our algorithm tries to achieve accurate boundary point detection by enhancing the boundary of interest. The core steps in our basic retinal layer boundary surface segmentation (BRLBSS) algorithm are below.
Step 1: Denoising A three-dimensional average filter with size of 1 2 3 K K K  is applied to the original retina OCT volume data to obtain a smoothed volume data, S , where 12 , KK and 3 K stand for the filter window width in the x, y, and z directions, respectively, and are all positive odd numbers.
Step 2: Extract the boundary points

1) Enhance the boundary with gradient information
A three-dimensional differential filter with a size of constructed using the pixel and its adjacent pixels, ,, Then, the differential filter is used to obtain a differential volume data, D. Specifically, for every cuboid in V, the intensity sum of the pixels above/below the cubiod center is subtracted from that of the pixels below/above the cuboid center, and the result is to obtain the intensity of the cuboid center.

2) Enhance the boundary by position, gradient and intensity
The depth, gradient and intensity information for the boundary are used to generate a boundary-enhanced volume data, I , to further highlight the boundary of interest. I is calculated as follows: where i, j, and k are the x, y, and z coordinates of each pixel, and 12 , ww are non-negative real numbers that are directly proportional to the depth coordinate, k .
The intensity of the desired boundary in I is likely to be the maximum value in its A-scan.

3) Extract the boundary surface points
In each A-scan of I , the point with the highest pixel intensity or the first peak (up to down or down to up, depending on the direction of the boundary of interest) is taken as the preliminary location of the desired boundary. Because the highest pixel intensity value is likely to lie on the desired boundary positions, the results obtained by this approach generally produce an accurate estimate of the location of the boundary.
Step 3: Correct the error points The depth positions with large errors are corrected using surface smoothness constraints, which require that the difference between the z coordinates of adjacent pixels is small. In The final depth information in A constitutes the desired boundary surface.

3． Implementation of the Algorithm for Segmenting Seven Retinal Layer Boundary Surfaces
This section details the implementation of the segmentation algorithm in Section 2 that automatically segments seven retinal boundary surfaces in SDOCT volume data.

RPE-Choroid boundary surface detection
For smoothing, the number of iterations also should be determined reasonably.
According    three-dimensional filters in the x direction (Fig. 2) should be used for the fovea region.

Vitreous-ILM boundary surface segmentation (VI_BSS) algorithm:
Step 1: Denoising For the flattened volume data V  , a three-dimensional average filter with window size 666  is used to produce a volume data denoted M . Then, a threshold filter is applied to M. In our experiment, the threshold value was 30. This filter procedure was iterated several times to obtain a smoothed volume data, S.
Step 2: Extract boundary points

1) Enhance the boundary with gradient information
A differential filter in the z direction is applied to the volume data S using a cuboid with a size of 1 1 11  . The new volume data is denoted by D.

2) Extract discrete boundary points
In every A-Scan of the volume data D, the first peak point is determined by moving from up to down. These points constitute the preliminary Vitreous-ILM boundary surface.
Step The VI_BSS algorithm is effective for OCT images with vessel shadows and/or a fovea region. However, it may fail for clinical OCT images with significant noise above the ILM (Fig. 9). To adapt this case, an improved algorithm is proposed.

Vitreous-ILM preliminary boundary surface segmentation (VI_PBSS) algorithm:
Step 1: Denoising Gray-scale morphological corrosion with a ball structure of radius r is performed on the smoothed volume data S obtained from the VI_BSS algorithm. Then, the average and threshold filter processing is repeated, as done in Step 1 of the VI_BSS algorithm. We denote the denoised volume data with SE. In our implementation, r = 5 was used.
Step 2: Extract the preliminary boundary points

1) Enhance the boundary with gradient information.
A differential filter in the z direction is applied to the volume data SE using a cuboid with a size of 1x1x11. The new volume data is denoted by D.

2) Extract discrete boundary points
In every A-Scan of the volume data D, determine the first peak point moving from up to down. To compensate for the excessive erosion, /2 r   pixels are subtracted from the depth coordinate, and these points constitute the preliminary Vitreous-ILM boundary surface points.
The VI_PBSS algorithm works well for OCT images with considerable noise above the ILM, but it may fail for clinical OCT images with prominent vessels (Fig. 10) due to excessive erosion.

ONL-IS/OS boundary surface segmentation
ONL-IS/OS boundary surface segmentation can be done with the flattened volume data V in between the RPE-Choroid and Vitreous-ILM boundary surfaces, where pixel intensities are set to 0 below the RPE-Choroid boundary and to 255 above the Vitreous-ILM boundary. We apply the three-dimensional differential filter in Eq. (2) with a size of 3 3 11  on the flattened volume data to obtain new volume data, D. In Eq. (3), 1 wk  and 2 0 w  are used to obtain an enhanced volume data, I. In every A-scan of data I, the point with the largest pixel value is taken as the preliminary location of the desired boundary. In principle, the ONL-IS/OS boundary surface can be obtained from step 3 in the BRLBSS algorithm. However, this algorithm may not work well for some OCT images with dark spots in the fovea region and near RPE, as shown in Fig. 12. z a a x a x a x     .
3) The x coordinate of every noise point is substituted into the polynomial z 1 , and the polynomial fitted value is taken as a depth estimate for the noise boundary point.
With the last obtained preliminary ONL-IS/OS boundary points, the error points are corrected using the same algorithm that is used for the RPE-Choroid smoothing.
The ONL-IS/OS boundary surface can be segmented correctly by using a polynomial fitting to eliminate the error points (Fig. 13).

OPL-ONL, NFL-GCL, IPL-INL, INL-OPL boundary surface segmentations
To segment the OPL-ONL boundary surface, we further shrink the search space in between the ONL-IS/OS and Vitreous-ILM layer boundaries. The basic steps are similar to the steps used in the ONL-IS/OS boundary segmentation. The differences are as follows: 1) The OPL-ONL boundary does not appear dark spots like the ONL-IS/OS boundary, and it may have a large curvature when the fovea region lies in the image. As a result, the step using a third-order polynomial to eliminate some error points is not used here.
2) The 3D differential filter in Eq.

4． Experimental results and analysis
To determine the segmentation accuracy of our algorithms for SDOCT volume data,  Table 1 shows that the average absolute error compared to manual delineation is 4.05μm, and the thickness of the thinnest layer in OCT volume is about 25 μm ( Table 2). Compared to the thickness, the absolute error is acceptable.  To further examine our algorithm, we calculated the thickness of the OCT layer from manual delineated data and auto-segmentation data. The experimental results are given in Table 2. The results listed in Table 1 and Table 2 are very similar in essential, which shows that our methods are rather practical for substituting the manual segmentation method. To further analyze the stability of our algorithm, we calculated the mean absolute errors for each OCT volume (Fig.15). It shows that the overall error for each subject is quite stable, but the segmentation of the layer INL-OPL is quite conditional on the quality of the data. Anyway, according to Table 1 and Table 2, these errors could be acceptable comparing to the thickness of the layer. OCT images with large curvature, severe noise above the ILM, fovea, prominent vessels, or large shadows because of diseases. The average computation time was 3.5 seconds per frame (Intel Core I7 CPU at 3.0GHz, and 8GB RAM), and the program can be optimized to reduce the time further. Therefore, our algorithm is simple and fast. Once the retinal layer segmentations have been performed successfully, every layer boundary surface can be visualized, and a thickness map between any two boundaries can be generated, which can be useful for disease diagnosis. An example visualization of the Vitreous-ILM and RPE-Choroid boundary surfaces is shown in Fig. 17, and the whole retinal thickness map determined from them is shown in Fig. 18.  Our method is inspired by the work in [22], but it also differs from it greatly. In principle, the method of Fabritius et al. only uses the pixel intensity of A-scans to detect the preliminary boundary points, and then boundary position information is used to eliminate the error points iteratively. As a result, their method is sensitive to noise and is invalid for images with severe noise. Our method uses pixel position information, gradient information and intensity information simultaneously to extract the boundary points so that the number of large error points is greatly reduced. Thus, the error point correction method used in our algorithm is completely different. Moreover, our algorithm can segment seven retinal boundary surfaces rather than just two. Compared to the complex 3D graph search-based methods in [18,19], our method is simple in principle and every efficient, and it does not depend on more prior knowledge and training dataset like the method in [20]. Additionally, our algorithm can deal with retinal OCT images with considerable noise above the ILM, and it does not need to detect prominent vessels or the fovea region explicitly. In our algorithm, the most three important boundary surfaces (i.e., RPE-Choroid, Vitreous-ILM and ONL-IS/OS) are segmented first, which can be used in improving some other algorithms such as the ones developed in [4,18,19] to shrink their search space at first, as those methods require a considerable amount of time to segment these boundaries because of the large search space.
As is the case for most of the current OCT image segmentation algorithms, our algorithm processes healthy or slightly abnormal retinal OCT images well but fails to accurately process retinal OCT images with serious diseases (Fig. 19).

5． Conclusions
We propose a novel 3D segmentation method for retinal OCT volume data that uses pixel intensity, boundary position information, and intensity changes on both sides of the layer borders simultaneously. The method designs a specific 3D differential operator for the processing boundary to enhance the border, it conducts a three-dimensional smoothing procedure to denoise the volume data, and it further utilizes the boundary position to produce an enhanced boundary volume data that serves as a better indicator for identifying the desired boundaries. Our method can segment seven boundary surfaces, and it is automatic, efficient and practical.