Automated retinal layer segmentation in optical coherence tomography images with intraretinal fluid

1. Introduction

The accurate segmentation of retinal layer boundaries in retinal images is a critical process for pathological analysis and diagnosis. However, the images of the retinal layers are sometimes complex and difficult to segment. For example, diabetes mellitus may lead to macular edema (ME).¹ In the optical coherence tomography (OCT) images of a retina affected by ME, a fluid-filled area can be seen; this area has low reflectance intensity and is a low contrast region in retinal tissue.² The fluid-filled area distorts the retinal tissue, which may cause a segmentation algorithm to fail.

So far, researchers have proposed a number of methods for retinal segmentation. Some methods are based on graph theory and dynamic programming.^3,4,5 Others are based on machine learning, such as pixel-wise classifiers.^6,7 Further, neural-network-based methods,^8,9,10,11 active contour methods, and level-set-based methods^12,13,14,15 have also been proposed. These methods exhibit good performance in segmenting images of healthy eyes, some of which have performed well even in cases involving retinopathy. Kugelman et al.’s⁸ and de Sisternes et al.’s⁹ methods can be used for eyes affected by age-related macular degeneration, whereas those of Vermeer et al.⁶ and Zang et al.¹¹ can be used to segment images of retinae affected by glaucoma. The methods of Lang et al.,⁷ Carass et al.,¹⁴ and Liu et al.¹⁶ have high accuracy in cases involving multiple sclerosis.¹⁷ These methods can handle some pathological eyes but not those affected by retinal fluid accumulation. Novosel et al.’s method¹⁵ can segment images of fluid-filled retinae affected by central serous chorioretinopathy.¹⁸ However, it is unclear whether the method can handle images of retinae affected by ME. In the eyes with retinal fluid accumulation, Wu et al. detected subretinal effusion,²⁷ Sappa et al.²⁸ and Hassan et al.²⁹ segmented and classified the retinal fluid accumulation. Wei and Peng,³⁰ Liu et al.,³¹ and Montuoro et al.³² segmented the retinal layer while separating the retinal fluid accumulation, achieving good results. But they rely on convolutional neural networks or unsupervised methods, which require a long time of training. Guo et al.’s method¹⁰ considered diabetic retinopathy, including ME, but the study did not cover the segmentation of a fluid-filled area because the structure of such an area is different from that of a healthy retina. In the dataset considered in this study, ME and microcystic ME sometimes appear simultaneously; hence, a segmentation method for images showing retinal fluid damage is needed.

Here, we propose a method to segment the retinal layer boundaries in OCT images with retinal fluid damage. Considering the structure of the retinal fluid region, fan filtering^19,20 and gamma correction were employed to enhance the linear information pertaining to the boundaries in retinal images. A random forest (RF) classifier.²¹ was used to predict the boundary probability for each pixel. The classifier used multiple decision trees to learn the retinal features and make predictions. In the final segmentation, boundary redirection (BR) and similarity correction (SC) were used to refine the predicted boundary. The proposed method segmented images of healthy and fluid-filled retinae found in our dataset into 11 boundaries (10 layers). The names (abbreviations) of these 11 boundaries are: (1) inner limiting membrane (ILM), (2) nerve fiber layer–ganglion cell layer (NFL–GCL), (3) ganglion cell layer–inner plexiform layer (GCL–IPL), (4) inner plexiform layer–inner nuclear layer (IPL–INL), (5) inner nuclear layer–outer plexiform layer (INL–OPL), (6) outer plexiform layer–outer nuclear layer (OPL–ONL), (7) external limiting membrane (ELM), (8) myoid zone–ellipsoid zone (MZ–EZ), (9) ellipsoid zone–outer segment layer (EZ–OSL), (10) outer segment layer–retinal pigment epithelium (OSL–RPE), and (11) Bruch’s membrane (BM).

2. Materials and Methods

2.1. Data acquisition

All the participants provided informed written consent. The study was conducted according to the Declaration of Helsinki. The protocol was reviewed and approved by the ethics committee. In this study, 52 OCT volume scans were acquired from one eye of 52 participants (including 40 healthy controls and 12 subjects with eyes containing intraretinal fluid; age range: 18–50 years) using a commercial OCT system (Topcon 3D OCT-1). Each volumetric data element consisted of 128 B-scans, each B-scan was composed of 512 A-lines, and each A-line included 885 sampling points. The volumetric scan captured retinal reflectance signals with a field-of-view size of $6\times 6\times 2.3$mm³ at the central macula region.

2.2. Preprocessing

Preprocessing consisted of two steps: filtering and flattening. In the filtering step, each B-scan was filtered by a fan filter with a passband angle $\theta =30^{\circ }$ to enhance the interlayer texture information.²⁰ Then, gamma correction was employed to enhance the image contrast. In the flattening step, which was the final process, the retinal boundaries were treated such that they were easier to find.

2.2.1. Fan filtering

The proposed approach used fan filtering, which can improve the interlayer texture information pertaining to the retinal image. A fan filter^19,20 is a frequency-domain filter, which was first used in seismic data processing. The fan filter has two common passband versions: horizontal and vertical passbands [Figs. 1(a) and 1(b)]. The gray area value is 1, and the white area value is 0. It can maintain directional information in a range by adjusting the passband angle $\theta$. As shown in Fig. 1(c), we tested the gradient enhancement effect of $\theta$ sector filter at 15–45^∘ (spacing of 5^∘), and its value first increased and then decreased. It has the maximum effect at 30^∘. The retina is mostly distributed horizontally. In order to improve the vertical contrast between different retinal layers, the retinal layer can be better segmented. A fan filter with a vertical passband was used in our method. In the frequency domain, the fan filter’s passband $H$ is defined by

where $u$ and $v$ are frequency variables and $\theta$ is the passband angle, as shown in Fig. 1(B). To confirm the effectiveness of the fan filter in retinal image processing, the frequency-domain images of OCT B-scans should be analyzed. For most B-scans in the spatial domain, the retinal layer is a curve-like band, and the layer boundaries can be distinguished by the contrast between layers. Linear information in the spatial domain corresponds to directional information in the frequency domain, which indicates that the linear texture of each layer in the spatial domain can be enhanced by reserving directional information [Figs. 1(d–i)].

2.2.2. Flattening

Flattening can shift most of the boundaries to horizontal, making them easier to be tracked. ILM and BM are used in flattening, as they have large boundary gradients (Fig. 4) and are easy to be found. Besides, the region of interest is between these two boundaries. In order to find the ILM and BM layers, we first performed fan filtering and gamma correction on each original B-scan to enhance the gradient and contrast of the original image and obtain the longitudinal gradient, as shown in Figs. 2(a) and 2(b). After Sobel filtering, it can be clearly found that the first three highest longitudinal gradient layers are ILM layer, MZ–EZ, and BM layer (Fig. 4). Therefore, the ILM layer is first found; the threshold interval is set according to the ILM boundary to shield the area. Then MZ–EZ is found, and finally, the BM layer is determined. Points that diverge from boundaries are removed due to noise and other factors. We use nearest-neighbor interpolation to complement these points on boundaries, and finally use Gaussian filter to smooth the ILM and BM boundary.

After obtaining the BM boundary, we calculate the offset distance between the BM layer and the central position of the image (at 885/2 pixels). Then the pixels between ILM and BM move the offset distance to the center position, which realizes retinal flattening at the central region of the image by the flattening operation, as shown in Figs. 2(c) and 2(d). It should be noted that the operations such as fan filtering and gamma correction made in Fig. 2 are only for better finding the ILM and BM layers and better flatness of the original image. Since the intensity of each B-scan is uniform, it may affect the segmentation results of the retinal boundaries. Therefore, the intensity of each B-scan should be normalized. The normalization method is max–min normalization, and the formula is as follows :

where $p_{\mathrm{\text{norm}}}$ is the normalized value; $p$ is the value of the original pixel; and max(B-scan) and min(B-scan) are the maximum and minimum values of the current B-scan.

2.3. Training and prediction

An RF classifier was trained to predict the retinal layer boundaries. To build the training dataset, 12 B-scans (six in the foveal region and six in the perifoveal region) were randomly selected from five subjects. The ground truth of the 11 retinal layer boundaries was manually labeled by experts.

In this study, 31 features were provided for RF classifier training⁷; these features helped the classifier to learn contextual, gradient, and position information pertaining to each class itself and the relationships between different classes. The details of the features are listed in Table 1. In particular, features 14–25 are generated by Gaussian anisotropy filtering. The filtering kernel is the full combination of three directions ($-10^{\circ }$, 0^∘, 10^∘) and two pixel scales $\sigma (x_1$, $y_1)=(5,1)$ and $\sigma (x_2,y_2)=(10,2)$. Features 12–25 are all made using the entire image.

In total, 31 features were used to train the RF classifier; a part of the features were set by referring to Lang et al.’s method.⁷ Because the retinal fluid cases had more complex structures than the healthy control cases, 27 features were sufficient for the healthy controls, whereas all the 31 features were required for the retinal fluid cases to achieve better accuracy. To identity healthy controls and subjects with eyes containing retinal fluid, a thickness threshold $T_m=370$$\mu$m was defined based on our retinal dataset. The neighborhood features for the healthy retina are shown in Fig. 3(a). Gradient features 12 and 13, which were used in the retinal fluid cases, were not used in the healthy control cases.

After training, the RF classifier can learn the features of the pixels from the training set. Accordingly, the RF classifier can predict the probability of the retinal layer boundaries in each pixel. Then, the retinal layer boundaries can be generated using probability maps [Fig. 3(b)].

2.4. Obtaining accurate retinal layer boundaries from probability maps

Each probability map records the probability of the corresponding boundary but not the exact position of each boundary. To track the boundary from the probability map, we used the combination of the BR and SC methods. These two methods are designed for different kinds of boundaries, and they exhibit good performance and strong robustness for the images of both healthy control and retinal fluid cases.

The proposed BR and SC methods both use nonmaximal suppression (NMS) and hysteresis thresholding in boundary tracking. NMS is used to determine the point with the maximum probability in each A-scan, and these points are divided into strong and weak boundary points. Hysteresis thresholding is used to keep the strong and weak boundary points attached to strong boundary points. When the boundary meets an intraretinal fluid region, the boundary probabilities are usually small, so the threshold for these kinds of boundaries also needs to be small.

The tracking of boundaries in a specific order is helpful to increase accuracy. In this study, the tracking order was chosen to be from the most reliable boundaries to the unreliable boundaries. The gradient curve of each boundary for retinal layers is shown in Fig. 4. It can be seen that the ILM layer, MZ–EZ, and BM layer have higher gradients, so they are considered as the reliable boundaries. First, these three boundaries are tracked, and then these three boundaries are used as restricted areas to delineate the other layers. The boundary tracking order is as follows: ILM, BM, MZ–EZ, EZ–OSL, OSL–RPE, ELM, NFL–GCL, GCL–IPL, INL–OPL, IPL–INL, and OPL–ONL.

For the BR and SC methods, the searching region settings in the same boundary were sometimes inconsistent because the thickness of a layer in each A-scan was not the same. Therefore, it was necessary to consider the thickness of each layer. We randomly chose 30 volumes from the dataset and segmented the retinal layers. Finally, the interval distance between the boundaries was calculated as the evidence for the searching region setting.

2.4.1. BR method

The BR method was used for tracking the boundary with a small curvature and a high prediction probability. This kind of boundary is common and easy to track, but problems still occur when mistaken predictions are carried out [Fig. 5(c)]. First, the BR method was used to obtain the preliminary boundary using NMS,^22,23 the probability map was then converted into a binary map through hysteresis threshold, and the cubic spline function was used to fit the binary map to obtain the fitted curve. Then, a new searching region was selected by choosing a region with a width of 5–10 pixels above and below the fitted curve in the first step, and the boundary was tracked again in the searching region. The result of the second round of tracking was considered as the final boundary. An example of the BR method is shown in Fig. 5. By using the BR method, the misjudgment [Fig. 5(c)] in the vessel shadow region can be corrected [Fig. 5(f)], as shown in the example.

2.4.2. SC method

The SC method is usually used for boundaries with a large curvature and a low prediction probability (e.g., NFL–GCL). These boundaries are difficult to track using the BR method. Thus, for these boundaries, boundary similarity (the shape of these boundaries is similar to the neighborhood boundaries) is the main factor to be considered in the searching region setting. For example, in B-scans near the fovea, NFL–GCL’s searching region is set by the ILM because it is the boundary above NFL–GCL. The ILM is the first tracked boundary in these 11 boundaries. Examples of using the SC method in the images of the retinal fluid cases are shown in Fig. 6.

3. Results

3.1. Algorithm performance

A large blood vessel in the superficial layer can produce hyper-reflection artifacts and strong shadow artifacts in the superficial and deep layers, respectively. These artifacts may cause the segmentation algorithm to fail; thus, we use fan filtering to eliminate these artifacts. Figure 7(b) shows the result of using a fan filter on the image in Fig. 7(a) with $\theta =30^{\circ }$. As shown in Fig. 7(b), the fan filter can suppress the hyper-reflection area and compensate for the shadow region (Fig. 7d). The width of the vessel shadow area can be well compensated by adjusting the passband angle of the fan filter. The larger the angle, the wider is the area from which the vessel shadow can be removed. In addition, gamma correction is used to enhance the contrast of the gradient on the layer boundaries.

For pathological cases, the fan filter can maintain the consistency of the boundaries in the lateral direction. The fan filter can also suppress the effect of pseudocysts and reduce the gradient level of the edema boundary in the lateral direction [Figs. 7(c) and 7(f)]. In general, the proposed filter can smooth the retinal image and reduce the gradient level in the lateral direction. The fan filter can further weaken partially discontinuous noise points and preserve the consistency of the layer boundaries. Moreover, the gradient on the retinal boundaries is enhanced by the fan filter.

The proposed boundary tracking method also shows good performance because the layer thickness provides more evidence about the searching region setting. In the large vessel shadow region, the segmentation usually fails [Fig. 7(g)] and the prediction is strongly influenced by the blood vessel. By adding the searching region, our algorithm can generate correct segmentation, as shown in Fig. 7(h). Furthermore, our method shows good results in boundaries affected by the retinal fluid region, as shown in Fig. 6.

3.2. Validation results

To validate the segmentation accuracy of the proposed method, we quantified the average errors between the manually delineated ground truth and auto-segmented results in the testing set, which included 40 healthy controls and 12 subjects with eyes containing intraretinal fluid. The comparison results are listed in Table 2. There is no statistically significant difference between the absolute errors of two groups, indicating that our classifier can accurately predict the location of retinal boundary of healthy control and subjects with diabetic macular edema (DME).

4. Discussion

In this paper, we proposed an auto-segmentation algorithm for retinal layer segmentation. The algorithm consisted of a novel filtering method to smooth the retinal volume and a new boundary tracking method. The proposed method can segment 10 retinal layers in the images of both healthy eyes and eyes with intraretinal fluid.

Filters such as the average filter and adaptive-weighted bilateral filter²⁴ were used for preprocessing. However, these filters cannot eliminate vessel shadow artifacts, which can cause the segmentation to fail. The fan filter can smooth the image by setting the passband angle, which corresponds to the linear features in a direction range. Compared with other filters, the fan filter is more effective in suppressing most vessel shadow artifacts. The average signed vertical gradient on the boundary obtained between the images filtered by the fan filter with different $\theta$ values and the original images [Fig. 2(c)] shows that in most cases, the boundary gradient in the filtered images is better than that in the original images. When the angle is set at 30^∘, the gradient of the layer boundaries in the fan-filtered image reaches a relatively high level.

In recent years, researchers have proposed numerous algorithms to segment retinal layers.^{3,8,10,16,25,26} The pixel classification method is more reliable than the gradient method because the former method considers not only the gradient but also other features such as the context and location. The gradient might be the most critical information in retinal layer segmentation, and some algorithms mainly rely on the gradient information to segment retinal layer boundaries.³ However, the role of the gradient is also important on the boundary of the fluid region, which may cause errors in segmentation. Recently, neural-network-based methods have become popular, and, according to several studies, their performance is better than those of other methods.^16,26 Compared with the neural network, the RF classifier is more controllable because the features can be selected manually. In addition, the RF classifier shows better performance than other kinds of classifiers,⁷ and hence, we chose the RF classifier in this study.

For boundary tracking, many algorithms have been used to handle probability maps, such as Lang et al.’s Canny and graph-based searching (GS) methods,⁷ Carass et al.’s multiple-object geometric deformable model method,¹⁴ and Liu et al.’s layer boundary evolution method.¹⁶ Compared to Lang et al.’s Canny method, their GS method has better accuracy but also requires greater computation time. Methods of Carass et al. and Liu et al. are both based on the level set method, which requires a higher computation capability for complex algorithms. Compared with other methods, the method proposed in this study has a lower computation complexity, so it needs a shorter running time. Additionally, it has a high accuracy because of searching region setting in the tracking stage. Thus, the proposed method may be a better solution for small training datasets. Although the searching region setting enables the algorithm to eliminate the influence of error prediction, it still requires prior knowledge about the dataset. Considering the heavy workload involved in manual segmentation, the thickness calculation required for the searching region setting can be realized through auto-segmentation and manual correction. The limitation of the proposed method is the lack of samples with a large retinal fluid size. Therefore, the performance of the proposed method may not be good for samples with a large amount of retinal fluid. In the future, the neural network will be considered to enhance segmentation accuracy, and the retinal fluid region will be segmented and analyzed to assist in diagnosis.

Conflict of Interest

The authors declare no conflict of interest.

Acknowledgments

This work was supported by Grants from the Research and Development Projects in Key Areas of Guangdong Province (2020B1111040001), the National Natural Science Foundation of China (NSFC) (81601534, 62075042, and 61805038), and Guangdong-Hong Kong-Macao Intelligent Micro-Nano Optoelectronic Technology Joint Laboratory (2020B1212030010).