1. Introduction

Glaucoma is a group of neurodegenerative eye diseases that are among the leading causes of irreversible blindness worldwide [1]. It is characterized by the degeneration of retinal ganglion cells (RGCs) and their axons in the peripapillary retinal nerve fiber layer (pRNFL). Additionally, it leads to specific alterations of the optic disc with loss of neuroretinal rim tissue, leading to slow progressive, irreversible loss of function. Although its pathogenesis is not completely understood, reduced ocular blood flow is suspected to contribute to its onset and progression [2–4]. If diagnosed and treated early, the progression can be slowed down or even stopped. Diagnosis of glaucoma is currently done by quantifying optic disc changes, visual field defects, and elevated intraocular pressure.

Fundus photography and scanning laser tomography have been one of the first image-based methods used to investigate glaucoma. The 2D en face images obtained from these devices allow us to examine the shape of the optic nerve head (ONH), for instance the cup-to-disk ratio or the rim area [5–9].

With the invention of optical coherence tomography (OCT) [10], three-dimensional in vivo imaging of the retina became feasible. This allowed for the visualization of structural information of the different retinal layers. Consequently, the thickness of the retinal nerve fiber layer (RNFL) [11,12] and the ganglion cell layer with inner plexiform layer (GCIPL) and ganglion cell complex (GCC = mRNFL + GCIPL) [13–15] were discovered as biomarkers for glaucoma. Moreover, in addition to the two-dimensional measurements of the ONH, three-dimensional measurements like its volume have been investigated [16,17].

With the development of OCT angiography, a relatively recent extension of OCT, a non-invasive, three-dimensional method for assessing retinal vasculature became available. This allows studying how glaucoma affects the microvasculature of the eyes. A widely used biomarker is vessel density (VD) which was studied in different regions (macula and ONH) and different sectors and plexus [18–24].

Since most forms of glaucoma are asymptomatic until the disease has progressed to an advanced stage, regular screening could detect the disease at an early stage. This would allow to intervene early and delay or even prevent the onset of symptoms. However, broad screening could result in an unmanageable workload for physicians. To prevent this, computer-aided diagnosis could be used to counteract this increase in measurements. Initial work on automatic screening used the previously identified biomarkers as computational features for machine learning to classify whether an eye exhibits signs of glaucoma. This was done by using RNFL and GCIPL thickness maps extracted from OCT volumes [25–27], microvascular features extracted from en face OCTA ONH scans [28], or a combination of VD and thickness parameters [29,30]. While for both OCT and fundus imaging fully data-driven approaches were presented [31–35], prior work using only OCTA still relies on handcrafted features.

To close this research gap we study the discriminative power of automatically detected features using convolutional neural networks (CNNs) for the task of distinguishing glaucoma patients from healthy controls and compare these features with handcrafted features from the literature. We evaluate the different approaches on 3x3 mm macular projections using the whole retina, the superficial vascular plexus (SVP), the intermediate capillary plexus (ICP), and the deep capillary plexus (DCP).

2. Methods

In this section, we will first explain our method of automatically extracting features using CNNs. After that, we introduce two sets of handcrafted features we used to benchmark our approach.

2.1 Automatically extracted features: convolutional neural networks

Different convolutional neural network (CNN) architectures were examined for automatic feature extraction. All of the models consist of a block of convolutional layers that extract the features from a two-dimensional input image followed by the actual classification. The networks used in this study are DenseNet121, DenseNet161 [36], ResNet-152 [37], ResNext-50-32x4d [38] and WideResNet-101-2 [39]. All of the models were pretrained on ImageNet. Since these pretrained models require the input to be a three-channel (RGB) image, the one-dimensional en face OCTA images were replicated three times such that the resulting images also had three channels. Moreover, the last layer was changed to have a single output indicating whether the patient suffers from glaucoma (label 1) or not (label 0). The training of these networks then consisted of two steps. In the first step, all layers, except the last one, were frozen and the network was trained until the validation loss did not decrease for three epochs. This was done because the last layer had to be randomly initialized in contrast to the frozen layers. By keeping the pretrained feature extraction constant, a classification based on these pretrained features was learned. Next, all layers were unfrozen and additional training was performed until a total of 200 epochs was reached. With this second step, the network could adapt the features and the classification to the different content of the input data. The loss was calculated using binary cross-entropy and optimized using stochastic gradient descent (SGD) with a momentum of 0.9. The batch size was set to 10, which was the maximum number of images that fitted into memory, and three different learning rates ($10^{-2}$, $10^{-3}$ and $10^{-4}$) were investigated. The mean and standard deviation of the pixel intensities of all images was computed before training and used to normalize the input images to have zero mean and a standard deviation equal to one. To deal with class imbalance, class weights were assigned to be inversely proportional to the total amount of images per class. This was necessary because our dataset had more en face OCTA images from glaucoma patients, which would bias the network towards classifying images as glaucomatous based on relative frequencies and thus lower the impact of actual pathological features in the image content. One drawback of CNNs is, that they can be considered a black-box approach. So it is be interesting to know which beneficial features were found on the 3 $\times$ 3 mm scans. To investigate this in greater detail, we applied the Gradient-weighted Class Activation Mapping (Grad-CAM) [40] algorithm to the predictions of the CNNs on the test set. The authors of Grad-CAM assume that the last convolutional layer of a CNN shows the best compromise between high-level semantics and detailed spatial information. They therefore use the gradient information flowing into this layer to assign importance values to each neuron for a particular classification target. These are then used to produce a class-discriminative localization map highlighting the important regions in the image for the prediction of a certain class.

2.2 Handcrafted features: Ong et al

To the best of our knowledge, there are few approaches known in literature for the distinction between healthy controls and glaucoma patients that use only OCTA data for this task. One example of handcrafted features that matched our dataset was presented by Ong et al. [28]. Although they presented their features for the region around the ONH, we adapted them and employed them on the macular region. In their work, they extracted the optic disc microvascular region for feature computation. Because the ONH is not included in our study, the features are directly applied to the macular region. They propose a set of local and global features for classification. The local features consist of Haralick’s information measures of correlation texture and features [41], inverse difference normalized and inverse difference moment normalized features [42], local structure mean, local structure standard deviation, and local structure deviation [28]. Global features included the mean, standard deviation, skewness, kurtosis, and entropy. They also proposed a thresholded cumulative count of microvascular pixels which was not used in this study because there was no microvascular region extracted.

After extracting these handcrafted features, we trained support vector machines (SVMs) to perform the classification task. Different kernels (linear, rbf and sigmoid) and regularization parameter (0.5, 1.0, 2.0 and 10.0) were tested. For the sigmoid kernel also different independent terms of the kernel function (0.0, 0.5, 1.0) were employed. As for the CNNs, class weights were used and the images were normalized in the same way.

2.3 Handcrafted features: vessel density

The vessel density (VD) is the most studied biomarker for glaucoma in OCTA images [18–21,24,43–47]. In order to compute the VD in this study, the histograms of all images were adjusted to match the histogram of a pre-selected scan in a first step. This procedure was done to neglect effects like different illumination across images such that the parameters of subsequent algorithms yield more reliable results once picked. After that, contrast limited adaptive histogram equalization (CLAHE) [48] is used in order to improve the contrast of the images, followed by the well-known Vesselness filter [49] to enhance the vessel structure in the images. As a last preprocessing step, hysteresis-thresholding was employed in order to create a binary image. This binary image is used to calculate the VD by dividing the number of pixels occupied by vessels by the total number of pixels in the image.

For the classification, we trained different SVMs on this feature with the same settings as described in the subsection above. This was done despite the fact that for a one-dimensional feature vector some of the parameters don’t matter and a SVM is unnecessary general. But since we also use multidimensional handcrafted features vectors for the evaluation, where a SVM is reasonable, we wanted to have the same learning strategy for the sake of an equal comparison.

3. Evaluation

In this section, we first describe the data used in this study and then the experiments performed on this data using the methods described in the section above.

3.1 Data

259 eyes of 199 subjects of the Erlanger Glaucoma Registry (Erlangen Glaucoma Registry, ISSN 2191-5008, CS-2011; NTC00494923) were analyzed retrospectively: 75 eyes of 74 healthy subjects and 184 eyes of 125 glaucoma patients. The cohort consisted of 98 male and 101 female subjects. The average age was 50.6 years with a standard deviation of 21.3 years. All subjects received a standardized ophthalmological examination including measurement of intraocular pressure (IOP) by Goldmann tonometry, fundus photography, and automated visual field testing. Measurements of retinal nerve fiber layer thickness (RNFL), retinal ganglion cell layer (RGC), inner nuclear layer (INL), and Bruch’s Membrane Opening-Minimum Rim Width (BMO-MRW) were done by Heidelberg OCT II Spectralis (version 1.9.10.0, Heidelberg Engineering, Heidelberg, Germany, Glaucoma Premium Module). En face OCTA imaging was done using Heidelberg Spectralis II OCT (Heidelberg, Germany). Images were recorded with a $15 {\times } 15^{\circ }$ angle and a lateral resolution of $5.7\;\mathrm{\mu}$m/pixel, resulting in a retinal section of $\sim 3 {\times} 3$ mm. The superficial vascular plexus (SVP), the intermediate capillary plexus (ICP), deep capillary plexus (DCP), and the whole retina (retina = SVP + ICP + DCP) were automatically segmented and projected using the manufacturer’s software. The projected enface OCTA images consist of 512 A-scans per B-scan and 512 consecutive B-scans.

The study has been approved by the ethics committee of the university of Erlangen-Nuremberg and performed in accordance to the tenets of the Declaration of Helsinki. Informed written consent was obtained from each participant.

3.2 Experiments

In this study, a five-fold cross-validation was performed. For one fold the data described above was split into 60% training set, 20% validation set, and 20% test set, with all eyes from one patient only belonging to a single set. The different architectures of the different methods as described in section Methods were then trained on the training set and evaluated on its validation set. The area under receiver operating characteristics (AUROC) was chosen as a metric for the validation. This procedure is then repeated five times such that each image appeared once in the validation set and once in the test set. After the five folds had been trained and evaluated, the mean AUROC per architecture over the validation results of the five folds was computed and the one, out of the three feature extracting methods, that yielded the highest mean AUROC was chosen for the final evaluation on the test set. This five-fold cross-validation was done for each of the three plexus and the whole retina.

4. Results and discussion

The best-performing model that was chosen for each method and plexus can be looked up in Table 1. The AUROC values on the test sets for the different methods and folds and the mean and standard deviation over all folds are displayed in Table 2.

Table 1. Architecture chosen per method and plexus (lr=learning rate, c=regularization parameter, coef0=independent kernel function term).

View Table | View all tables in this article

Table 2. AUROC values on the test sets for the different methods and folds and the mean and standard deviation over all folds.

View Table | View all tables in this article

Here we can see several things. First of all, the CNNs and the features proposed by Ong et al. performed better than the VD across all plexus. This might be because 3 $\times$ 3 mm scans were used and several studies showed that the most vulnerable areas for glaucoma are mostly outside the 3 $\times$ 3 mm area but inside the 6 $\times$ 6 mm area and thus the larger field sizes yield higher AUROC values than the small ones [46,47]. Moreover, on our data, the VD on the ICP performed best for the discrimination between glaucoma patients and healthy controls. Takusagawa et al. [18] reported on 6 $\times$ 6 mm scans an AUROC of 0.96 for the VD computed from the SVP. They also reported that on their data the VD was only slightly decreased for the ICP and not decreased on the DCP. Rabiolo et al. [50] evaluated different methods for the binarization of the vessel structure on 6 $\times$ 6 mm and reported an AUROC value of 0.94 on the SVP and of 0.99 on the DCP. But there were no exact values for the VD in the ICP found in the literature.

The CNNs and the SVM trained on the features presented by Ong et al. performed best on the SVP. Moreover, we can see that both yield approximately the same results for the ICP (CNN: 0.915; SVM(Ong): 0.912) and the DCP (CNN: 0.910; SVM(Ong): 0.917) but the CNNs outperform the SVM(Ong) approach on the SVP (CNN: 0.967; SVM(Ong): 0.942) and especially on the whole retina projection (CNN: 0.923; SVM(Ong): 0.789).

An interesting point here is that the SVM(Ong) approach yielded weak results on the whole retina projection in comparison to its mean AUROC values for the different plexus. A reason for this behavior might be because in the projections of the different plexus, individual vessels are better distinguishable from each other and are not superimposed over each other. Since the features proposed by Ong et al. are mainly textural features, this clearer visualization of the vessel network might be beneficial for automatic classification.

Table 2 shows that the CNN using the SVP achieves the highest AUROC values, which are comparable to the values reported in the literature for larger and/or other regions of the human retina.

Representative visual examples for the Grad-CAM heatmaps for the classification of the glaucoma class for correctly classified samples and falsely classified samples can be found in Fig. 1 and Fig. 2 respectively. There has been some debate on how reliable the interpretation of these saliency maps are [51]. However, when looking at Fig. 1 and Fig. 2 it seems plausible that in our case they highlight the regions in the images that contribute to classifying an image as glaucomatous. First of all, in the images that were classified as not exhibiting signs of glaucoma (e.g. being healthy) there are no or only small regions that are highlighted. On the other side, if an image was classified as glaucomatous, the number of highlighted regions increased and are larger. Assuming that these highlights show the regions that influence the classification decision, it is striking that the highlights seem to be in regions where vessels are sparser. Moreover, the network seems to focus mostly on the corners and borders of the images, which is in accordance with the aforementioned studies that showed that the vessel density works better in regions farther from the macula. This leads us to the conclusion that vessel density as a feature does not perform worse in general on the smaller scan sizes, but rather the regions in which it was calculated on were too large and/or did not fit. Consequently, improving the region selection for the computation of the vessel density might improve its performance for classification. Further research is warranted to investigate whether these conclusions drawn from the saliency maps hold.

 

Fig. 1. Grad-CAM heatmaps for correctly classified samples in the macular area of the test set. Darker/brighter colors indicate regions with lower/higher attention. No color indicates that there was no attention calculated by the Grad-CAM algorithm.

Download Full Size | PDF

 

Fig. 2. Grad-CAM heatmaps for falsely classified samples in the macular area of the test set. Darker/brighter colors indicate regions with lower/higher attention. No color indicates that there was no attention calculated by the Grad-CAM algorithm.

Download Full Size | PDF

There are several limitations to our study. We did not differentiate between different sub-types and severity stages of glaucoma. So it would be interesting to apply the automatic feature extraction to these more challenging tasks. Another point to notice is that we only discriminated glaucoma patients from a healthy control group and did not include other diseases like age-related macula degeneration (AMD) or diabetic retinopathy (DR) in this second group. It would be an interesting study to see how the different features perform in this more realistic setup. However, the Erlanger glaucoma registry that we worked with only provided data for glaucoma patients and the healthy control group but no other diseases. Another limitation is that for the vessel density, several studies have shown that the 6 $\times$ 6 mm region is more susceptible to changes. Our study was performed on the smaller 3 $\times$ 3 mm images because a larger amount of scans was available. Therefore, it would be interesting whether the features proposed by Ong et al. and the automated features extracted from the CNN could yield even better results if trained on the larger 6 $\times$ 6 mm regions, or if this is only beneficial for classification based on vessel density. However, we have seen that the AUROC values using the CNNs on the smaller scan size are already comparable to the ones achieved by the vessel density on the larger scan size. Consequently, it is feasible to use 3 $\times$ 3 mm scan sizes for the glaucoma detection on OCTA images when choosing appropriate features. Since CNNs are a black-box approach we can conclude that the automatically derived features in the 3 $\times$ 3 mm area allow a good distinction between glaucoma patients and healthy controls, but those features are difficult to interpret. So it would be interesting as a future work to investigate, especially in the SVP, what changes appear and if those could be used as an interpretable biomarker. We did a first step in this direction by computing saliency maps on the test set, but further research has to be conducted in order to demonstrate that the conclusions drawn from them are correct.

Overall this study shows that the automatically extracted features perform no less and in certain plexus even much better than the handcrafted features available for OCTA images for the distinction between glaucoma patients and healthy subjects on 3 $\times$ 3 mm macular scans. With these automatically extracted features we achieved AUROC values that are comparable to the best values reported in the literature, although using the smaller 3 $\times$ 3 mm scan sizes that have performed considerably worse in previously published studies.

5. Summary and conclusion

In this study, we examined handcrafted and automatically trained features for the classification of glaucoma patients. We show that the best performance was achieved on the SVP using the automatically extracted features. Moreover, AUROC values can be achieved that are comparable to the best-reported values in the literature, even on 3 $\times$ 3 mm scan sizes, which were shown to perform worse than larger scans for the vessel density, because the area of the macula that is most vulnerable to glaucomatous vascular damage is outside the central 3 $\times$ 3 mm region. Further analysis of the CNN using attention maps showed that the network focuses on areas where there are fewer vessels, but that these regions are mostly at the border and the corners of the scans and therefore can not be picked up so well by the mean of the vessel density over a larger area.