Automatic Three-dimensional Detection of Photoreceptor Ellipsoid Zone Disruption Caused by Trauma in the OCT

Detection and assessment of the integrity of the photoreceptor ellipsoid zone (EZ) are important because it is critical for visual acuity in retina trauma and other diseases. We have proposed and validated a framework that can automatically analyse the 3D integrity of the EZ in optical coherence tomography (OCT) images. The images are first filtered and automatically segmented into 10 layers, of which EZ is located in the 7th layer. For each voxel of the EZ, 57 features are extracted and a principle component analysis is performed to optimize the features. An Adaboost classifier is trained to classify each voxel of the EZ as disrupted or non-disrupted. Finally, blood vessel silhouettes and isolated points are excluded. To demonstrate its effectiveness, the proposed framework was tested on 15 eyes with retinal trauma and 15 normal eyes. For the eyes with retinal trauma, the sensitivity (SEN) was 85.69% ± 9.59%, the specificity (SPE) was 85.91% ± 5.48%, and the balanced accuracy rate (BAR) was 85.80% ± 6.16%. For the normal eyes, the SPE was 99.03% ± 0.73%, and the SEN and BAR levels were not relevant. Our framework has the potential to become a useful tool for studying retina trauma and other conditions involving EZ integrity.

Scientific RepoRts | 6:25433 | DOI: 10.1038/srep25433 on only a single 2D cross section image. In ref. 10 and ref. 18, the partial OCT projection image, or en face image, was developed to better visualize a map of the photoreceptor integrity and its disruption. However, the measurement of the EZ disruption area was still based on a 2D image and a manual method, which could involve subjective factors when selecting the disruption margins. Sayanagi et al. 11 developed an automated EZ disruption margin detection and a weighted EZ disruption area calculation method. However, their measurement was based on the assumption that the disruption region was circular, while in fact this region may have an irregular shape. Itoh et al. 17 developed an automated EZ mapping tool to assess the EZ integrity by providing en face visualization of EZ integrity, EZ-retinal pigment epithelium (EZ-RPE) height and EZ-RPE volume.
In this paper, we propose an automatic 3D framework to detect EZ disruption in macular SD-OCT scans. Machine learning classifiers have been widely used in a variety of OCT specific applications including layer segmentation 19,20 , drusen classification 21 , and microcystic macular edema segmentation 22 . In this paper, we apply an adaptive boosting (Adaboost) [23][24][25] -based method to classify the pixels as disrupted or non-disrupted. The contributions of this work are as follows: (1) we propose a novel framework for a 3D, automated, and quantitative analysis of EZ integrity in retinal OCT images; and (2) because the detection of EZ disruption is a typical imbalanced classification problem 26 , we apply two strategies on two different levels to improve the classification performance. These two strategies include the Adaboost ensemble classification algorithm at the algorithm level and an under-sampling strategy at the data level.

Results
Data Analysis. To evaluate the performance of the proposed method, all the EZ disruption regions in the 3D SD-OCT images were manually marked by an ophthalmologist slice by slice using the ITK-SNAP software 27 and saved as the ground truth. The leave-one-out method was used to train the Adaboost based integrated classifier models. Because the sample ratio of the majority class (non-disrupted) and the minority class (disrupted) was approximately (110 ± 256):1 on average, non-disrupted samples were randomly selected to match the disrupted ones. The EZ disruption volume was calculated by multiplying the disruption number by the voxel resolution.
The mean and 95% confidence intervals of the segmented EZ disruption region volumes were compared between eyes with retinal trauma and normal eyes. Student's t-test was used to evaluate the statistical significance of the disruption volume differences between the two groups of eyes. Statistical correlation analysis and Bland-Altman plot analysis were utilized for a performance comparison between the proposed method and the ground truth.
To assess our experiments, several measures based on the segmented volume of the EZ disruption including sensitivity (SEN), specificity (SPE) and balanced accuracy rate (BAR) were adopted. These evaluation indexes are commonly used in imbalanced classification problems and are defined as below: where TP, FN, TN and FP represent true positive, false negative, true negative and false negative, respectively.
Experimental Results. Figure 1 shows one of the detection results (Case #4 in Table 1) using the proposed framework, and the corresponding ground truth for the EZ disruption region. The en face projections of the original VOIs, ground truth, and corresponding detected EZ disruption are also shown in Fig. 1 Fig. 3. Student's t-test demonstrated a strong statistical significance for the detected EZ disruption volume differences between the two groups of eyes (p = 9.9112 × 10 −8 ≪ 0.001). Table 1 shows the detected EZ disruption volume, ground truth volume, whole EZ volume, SEN, SPE and BAR for the 15 eyes with retinal trauma. For the eyes with retinal trauma, the SEN was 85.69% ± 9.59%, the SPE was 85.91% ± 5.48%, and the BAR was 85.80% ± 6.16%. For the normal eyes, the SPE was 99.03% ± 0.73%. Because there were no true positives, the values of SEN and BAR were irrelevant.
For the eyes with retinal trauma, the correlation between the segmented EZ disruption volume and the ground truth was r = 0.8795 with a significance level p < 0.0001. The 95% confidence interval for r was 0.6683 to 0.9595. Figure 4 shows the Bland-Altman plot for the consistency analysis between the automatic segmented EZ disruption volume and the ground truth. We can see from Fig. 4 that they are consistent.

Discussion and Conclusion
In this study, we developed and evaluated an automatic method to detect the 3D integrity of the EZ in eyes with retinal trauma. Because the disrupted voxels in the EZ region are much less numerous than the non-disrupted Scientific RepoRts | 6:25433 | DOI: 10.1038/srep25433 ones, this leads to a typical imbalanced classification problem. To overcome this problem, an Adaboost algorithm (at the algorithm level) and dataset balance strategies (at the data level) are utilized. The vessel silhouettes and isolated points are excluded to decrease the false detections, using a vessel detector 28 and morphological opening operations, respectively. The average detected EZ disruption volume in the eyes with retinal trauma was statistically much larger than the corresponding volume in the normal eyes (Student's t-test, p = 9.9112 × 10 −8 < 0.001). In the eyes with retinal trauma, the SEN, SPE and BAR, using the proposed method, were 85.69% ± 9.59%, 85.91% ± 5.48%, and 85.80% ± 6.16%, respectively. There was a strong correlation between the segmented EZ disruption volume and the ground truth (r = 0.8795). In the normal eyes, the SPE was 99.03% ± 0.73%. This study has several limitations. (1) Although many studies have shown that the disruption extent of the EZ is an important clinical indicator for the injury degree of the photoreceptors and that EZ disruption may be closely associated with visual acuity in different eye diseases [7][8][9][10][11][12][13][14][15][16][17] , there are some controversies [29][30][31][32] . Whether the quantitative disruptions of the EZ have quantitative relationships with visual acuity and visual outcome has not yet been addressed. This study focuses on quantitative measurements of the EZ disruption volume; a further study will be carried out in the near future to determine the quantitative correlation, if any, of these quantitative measures to visual acuity and to the outcome of eyes with retinal trauma. (2) The sensitivity and specificity of the   proposed method could also be further improved. The classification errors could be due to two reasons. i) The inaccuracy of the surface segmentation results during the image pre-processing stage may cause unreasonable extractions of VOIs. For example, Fig. 5a shows an original B-scan, and Fig. 5b shows its segmented 7 th and 8 th surfaces; both of which have lower values than the correct surfaces on the left side. Figure 5c shows the ground truth (in red), and Fig. 5d shows the detected results (in yellow). It is obvious that there are both false positives and false negatives on the left side. ii) Because of the poor quality of the SD-OCT images, the ground truths marked by the ophthalmologist are subjective, especially in the transitional region between the disrupted and the non-disrupted regions. For example, in Fig. 5e, the quality of the B-scan is very low from the left side to the central fovea area. It is difficult to decide whether the corresponding EZ is disrupted or not. Figure 5f shows the ground truth (in red), which may contain subjective preferences. Figure 5g shows the detected EZ disruption results (in yellow), which are not very consistent with the ground truth. (3) Due to the collection difficulties and poor quality of the image data, in this paper, the proposed method was tested on a sample size of 30 images dataset (15 eyes with retinal trauma and 15 normal eyes). The testing dataset is not large, however, we are still collecting more data and will validate our method on a larger dataset in the near future. (4) This work is focused on the 3D shape and volume of the EZ disruption in SD-OCT image, for more comprehensive EZ disruption analysis, the retinal information from other imaging modalities such as fundus and fundus fluorescein angiography image should be considered as well.

Study Subjects and Data Collection. The Institutional Review Board of the Joint Shantou International
Eye Center approved this study and waived informed consent due to the retrospective nature of this study. Our study also complies with the Declaration of Helsinki. The patient records/information was made anonymous prior to analysis. The medical records and OCT database of JSIEC were searched and reviewed from February 2009 to December 2012. In total, 15 eyes in subjects with retinal trauma and 15 eyes in normal subjects were included and underwent a macular-centred (6 × 6 mm) SD-OCT scan (Topcon 3D OCT-1000, 512 × 64 × 480 voxels, 11.72 × 93.75 × 3.50 μ m 3 , or 512 × 128 × 480 voxels, 11.72 × 46.88 × 3.50 μ m 3 ). There were 12 males and 3  females in the trauma group, with a mean age of 30.3 ± 11.3 years (range: 8-43 years). There were 9 males and 6 females in the normal group, with a mean age of 33.1 ± 10.8 years (range: 7-46 years). Subjects with other eye diseases were excluded except for those with refractive error < = + /− 6 diopter. The raw, uncompressed data were exported from the OCT machine in ".fds" format for analysis.

Overview of the EZ Disruption Detection Method. The proposed method consists of three main
parts: pre-processing, classification, and post-processing. In the pre-processing step, the SD-OCT images are first denoised using fast bilateral filtering 33 . Then, the SD-OCT volume is automatically segmented into 10 intra-retinal layers using the multi-scale 3D graph-search approach [34][35][36][37][38][39][40] , which produces 11 surfaces. The retina in the original SD-OCT volume is flattened, and the 11 th surface (the bottom of the retinal pigment epithelium) is used as the reference plane. The EZ region between the 7 th and 8 th surfaces is extracted, which is the volume of interest (VOI) for our analysis. In the classification step, five categories, from a total of 57 features, are extracted for each voxel in the VOIs. Then, principle component analysis (PCA) is adopted for feature selection. Because the disrupted voxels (the minority) in the VOIs are far less numerous than the non-disrupted ones (the majority), it is a typical imbalanced classification problem. To improve the performance of the classification, we apply the following two strategies in the classification training: (1) an Adaboost algorithm is adopted to train some weak classifiers into an integrated strong classifier at the algorithm level; and (2) the majority samples are randomly under-sampled at the data level. In the classifier testing step, every voxel in the VOIs is classified as disrupted or not disrupted. In the post-processing step, the blood vessel silhouettes are identified and excluded by a vessel detector 28 and the isolated points are excluded by morphological operations to avoid false detections. Finally, the volume of the disrupted EZ is calculated. Denoising by Bilateral Filtering. Speckle noise is the main noise in OCT images 41 , and it affects the performance of image processing and classification. In this paper, we propose applying the bilateral filtering 42 method for denoising because it can remove speckle noise from images effectively while maintaining edge-like features. We have used a fast approximation algorithm 33 to reduce the computation time without significantly impacting the bilateral filtering result. Each B-scan (X-Z image) of the OCT images is smoothed separately by bilateral filtering.

Segmentation of the Intra-retinal Layer and the EZ Region. The filtered SD-OCT volume is then
automatically segmented into 10 intra-retinal layers using the multi-scale 3D graph-search approach [34][35][36][37][38][39][40] , which produces 11 surfaces (see Fig. 6). Then, all the surfaces are smoothed using thin plate splines. The retina in the original SD-OCT volume is flattened by adjusting the A-scans up and down in the z-direction, where the 11 th surface (the bottom of the retinal pigment epithelium) is used as a reference plane because of its robustness. Then, the EZ regions between the 7 th and 8 th surfaces are extracted as the volumes of interest (VOIs).
Feature Extraction. For classification, the following five types of low-level features are extracted: normalized intensity, block mean, block standard deviation, the absolute intensity difference in the 13 directions to be described later (step = 1, 2), and the grey-level co-occurrence matrix (GLCM) based features (contrast, correlation, energy and homogeneity in 13 directions). Therefore, 57 features are extracted, which are listed in Table 2.
The normalized intensity represents the voxel's grey level. As shown in Fig. 6, the intensity level of the disrupted region of the EZ is lower than the intensity level of the non-disrupted region. Therefore, if a voxel's normalized intensity level is low, it has a higher probability of being classified as a disruption, and vice versa.
The block mean and block standard deviation represent the average intensity level and the variance of the intensity level, respectively, in the local region centred around the voxel (region of 5 × 5 × 5 voxels). The absolute intensity difference in 13 directions represents the variance of the intensity between the centre voxel and its neighbours in 13 directions. Let α 1 stand for the angle between the X-axis and the projection direction on the X-Y plane, and let α 2 stand for the angle between the X-Y projection direction and the Z-axis. The 13 directions can be described as follows: (α 1 ,  The grey-level co-occurrence matrices (GLCMs) for the 3D volumetric data describe the spatial dependence of grey levels across multiple slices 43,44 . The 3D method searches for other grey levels in the 13 directions (mentioned above) across multiple planes and constructs 13 GLCMs. Here, the GLCMs in the 13 directions of every 5 × 5 × 5 block are constructed. Then, the following four features are calculated: (1) the contrast, which measures the local contrast of the volumetric image and is expected to be higher when a large grey-level difference occurs more frequently; (2) the correlation, which provides a correlation between the two voxels in a voxel pair and is expected to be higher when the grey levels of a voxel pair are more correlated; (3) the energy, which measures the number of repeated voxel pairs and is expected to be higher if the occurrence of repeated voxel pairs is higher; and (4) the homogeneity, which measures the local homogeneity of a voxel pair and will be larger when the grey levels of each voxel pair are more similar. Feature Selection. Based on above definitions, we have a total of 57 features extracted for each voxel in the VOIs. To reduce the dimensionality of the feature vector and describe the inter-correlated quantitative dependence of the features, a feature selection procedure based on the PCA is performed. In our experiments, the first 10 principle components are selected as the new features; they represent more than 90% of the information in the original features.
Adaboost Algorithm and the Under-sampling Based Integrated Classifier. In this study, the number of non-disrupted samples in the EZ region is far greater than the number of disrupted ones. The disrupted EZ samples and non-disrupted EZ samples belong to the minority and majority classes, respectively. This is a typical imbalanced classification problem, which means the class distribution is highly skewed. Most traditional single classifiers, such as the support vector machine, the k-nearest neighbour classifier, quadratic discriminate analysis, and the decision tree classifier, tend to show a strong bias towards the majority class and do not work well for this type of problem because they aim to maximize the overall accuracy. The Adaboost algorithm based integrated classifier [23][24][25] is one solution to overcome this problem at the algorithm level; it integrates multiple weak classifiers into a strong classifier and is therefore more sensitive to the minority. Hence, the Adaboost algorithm is adopted in this study.
To further improve the classification performance at the data level, the training datasets are balanced by under-sampling majority samples. In the training step, the Adaboost algorithm-based classifier model is calculated according to leave-one-out cross-validation, using all the disrupted samples and an equivalent number of randomly selected non-disrupted samples. In the testing stage, each voxel in the VOIs is classified as disrupted or non-disrupted using the trained Adaboost model.  voxels between the EZ and RPE are selected and each pixel in the 2D projection image is the average in the z-axis direction of the selected voxels at that particular x, y location in the OCT volume. Then, the vessel silhouettes are segmented using a KNN classifier. If the detected EZ disruption regions have the same x and y location as the vessel silhouettes, these regions are regarded as normal and removed as false detections. Due to the physiological connectivity of the EZ disrupted/non-disrupted regions, isolated disrupted/non-disrupted voxels are eliminated through morphological opening operations, where the shape of the structural element is set as ball with a radius of 5 voxels.