Unsupervised feature extraction of anterior chamber OCT images for ordering and classification

We propose an image processing method for ordering anterior chamber optical coherence tomography (OCT) images in a fully unsupervised manner. The method consists of three steps: Firstly we preprocess the images (filtering the noise, aligning and normalizing the resolution); secondly, a distance measure between images is computed for every pair of images; thirdly we apply a machine learning algorithm that exploits the distance measure to order the images in a two-dimensional plane. The method is applied to a large (~1000) database of anterior chamber OCT images of healthy subjects and patients with angle-closure and the resulting unsupervised ordering and classification is validated by two ophthalmologists.

Unsupervised feature extraction of anterior chamber oCt images for ordering and classification pablo Amil 1 , Laura González 2 , elena Arrondo 2 , Cecilia salinas 2 , J. L. Guell 2 , Cristina Masoller 1 & Ulrich parlitz 3 We propose an image processing method for ordering anterior chamber optical coherence tomography (oCt) images in a fully unsupervised manner. the method consists of three steps: Firstly we preprocess the images (filtering the noise, aligning and normalizing the resolution); secondly, a distance measure between images is computed for every pair of images; thirdly we apply a machine learning algorithm that exploits the distance measure to order the images in a two-dimensional plane. the method is applied to a large (~1000) database of anterior chamber OCT images of healthy subjects and patients with angle-closure and the resulting unsupervised ordering and classification is validated by two ophthalmologists.
Machine learning methods are extremely useful in biomedicine 1,2 and in particular for glaucoma detection [3][4][5][6] . The use of such methods can help to optimize the available human resources, to increase accuracy in diagnosis, and to make treatment decisions faster.
Glaucoma is the leading cause of global irreversible blindness 7 . Early diagnosis and treatment is a challenge given that glaucoma presents no symptoms in its early stages 8 . Diagnostic of angle-closure is based on the clinical observation of the angle at the slit-lamp requiring a goniolens that is placed on the patient's cornea. Anterior Segment Optical Coherence Tomography (AS-OCT) is a fast, useful and contact-less tool that allows visualization and measurement of the anterior chamber angle [9][10][11] . Various techniques exist and are being developed to improve image quality [12][13][14][15] , and work is also focused on the development of advanced tools to analyze such images [16][17][18] . However the quality of the images is not always enough for an accurate diagnosis and for example 19 , two glaucoma experts working together were unable to locate the scleral spur (see Fig. 8) in 28% of the images. Such landmark is of utmost importance to determine most of the relevant features that are used in angle-closure diagnosis.
Previous approaches for processing anterior chamber OCT images used segmentation algorithms and further analysis, while others used manual landmark determination. In works by Tian et al. 16 and Wu et al. 20 segmentation algorithms for OCT anterior chamber images were proposed, both discuss the difficulties of such task due to the noise in the image and to other artifacts. In 21 the authors analyzed data from manual landmark determination to classify images into five subcategories of angle-closure glaucoma. They used both supervised and unsupervised feature selection and AdaBoost classifiers (supervised learning) to achieve accuracies in the range of 84~87%. In 22 a method was proposed to extract almost 3000 features from the raw images, then, the most relevant features for classification were supervisedly selected, and machine learning was applied to the selected features. However, the fact that only 74 images were analyzed while 15 features were used, compromised the statistical significance of the study. In 23 Histogram of Oriented Gradients (HOG) features of anterior chamber OCT images and Support Vector Machine (SVM, supervised learning) were used to classify different glaucoma subtypes, achieving accuracies in the order of 80%.
In this paper we propose an image processing method that unsupervisedly orders OCT anterior chamber images according to what we demonstrate to be relevant extracted features. The reliability of the algorithm is tested with a large number of images (~1000). Importantly, our method is fully autonomous and can be used to analyze images with a wide spectrum of quality, even those with high levels of noise and artifacts.

Results
The outcome of the algorithm applied to the image database (both described in Sec. Methods) is presented in the way of an "Image Map". A regular grid in the coordinates space (w, v) is defined and one image per grid point is displayed. The results are presented in Figs 1-4 that show the Image Map obtained after applying IsoMap to the Euclidean distance ( Fig. 1), to the aligned Euclidean distance (Eq. 4, Fig. 2), to the aligned Hellinger distance (Eq. 3, Fig. 3) and after applying t-SNE to the Hellinger distance (Fig. 4). These figures are ordered with increasing complexity and performance of the corresponding algorithm.
In Fig. 1 (IsoMap with Euclidean distance) it is apparent that the algorithm ordered the images according to the orientation (horizontal axis) and position (vertical axis) of each eye inside the OCT image, which are irrelevant features, this simple algorithm is, then, not capable of extracting useful features. However, when including the alignment step in the pre-processing of the images which removed the variability detected by the first algorithm, meaningful features were extracted, as it can be seen in the Image Map in Fig. 2 (these features correlate with the features derived manually, as shown in Table 1). A similar map was obtained with the Hellinger distance, shown in Fig. 3 which marginally improves the performance with the Euclidean distance. To test the robustness of the image ordering, the t-SNE algorithm was applied (instead of IsoMap). The map obtained is shown in Fig. 4, which turned out to perform slightly better than IsoMap according to the correlations shown in Table 1.
In Fig. 5 the features returned by the unsupervised algorithm, using t-SNE, are compared with the features obtained from the manual annotation of the images: in the left panel the color code indicates the chamber depth, and in the right panel, it indicates the mean angle (average of α and β) of each annotated image. Clearly, the features obtained from the manual annotation correlate very well with the features returned by the unsupervised algorithm. However, one should notice that the annotated features are not independent, but strongly correlated with each other. A summary of the correlation between the ordering (Mappings) and different features is shown in Table 1.
Finally, the manual classification done by two expert ophthalmologists is included in the t-SNE (w, v) map. In Fig. 6, left panel, the color code indicates the different classes, and it can be noticed that the data points that represent images in the "wide-open" category are scattered in the left-down corner of the map, while the data points that represent images in the "closed" category are scattered in the right-top corner of the map.
In order to demonstrate that the separation between different classes obtained from features retrieved from the manual annotation of the images is comparable to the one obtained from the t-SNE features, the manual annotation features (mean angle and smallest ARA_500) are plotted in Fig. 6, right panel. A comparison with the left panel reveals that the manual annotation features also do not allow for a clear-cut separation of the different classes, and therefore, for classification purposes, the features obtained with the t-SNE map are as good as the manual annotation features.

Discussion
We have proposed a new algorithm for ordering anterior chamber OCT images in such a way that it is possible to classify them, in a fully unsupervised manner, in meaningful groups according to relevant features. We have tested the algorithm with a large set of images classified by two expert ophthalmologists, and with a larger set of annotated images. We have verified that the separation in the different classes defined by the ophthalmologists (closed, narrow, open, and wide open) is similar when using the manually extracted features, or when using the features that are returned by the unsupervised algorithm (Fig. 6).
Therefore, the abstract features generated by the algorithm provide novel tools for assessing OCT images of the anterior chamber. They can be used for direct classification of the images and, furthermore, they can be linked to established quantities used for characterizing diseased eyes (like chamber depth, iris-corneal angle) resulting in an automatic detection system. As the algorithm is fully unsupervised, it can be easily automated and set up in OCT imaging systems to aid technicians and doctors in an early diagnosis.
The two main advantages of the algorithm demonstrated here over previous works are that it doesn't need any ground truth or gold standard for training, and it does not rely on specific landmarks; thus, it can analyze images in which relevant landmarks are not visible or not easy to locate.

Methods
Data. The data consists of 1213 OCT images taken from consenting patients at IMO (Ocular Microsurgery Institute, Barcelona) 24 . The images were acquired using a Visante OCT instrument (Carl Zeiss Meditec) using the "Anterior segment" scan and the "Enhanced anterior segment" scan. The original resolution of all the images is 256-by-1024 pixels corresponding to an area of 16 mm-by-8 mm. All the patients included in the study consented that their data and images could be used for research and teaching purposes and the study has been approved by the ethical committee for clinic research in IMO. The patients were selected based on the already available data at IMO (retrospective study). We selected 247 patients at IMO database, 81 of which were glaucoma patents and the rest was a mixed group of healthy patients and patients with other diseases different from glaucoma and with normal intraocular pressure, 104 of those were cataract patients, from which we used images from before and after intraocular lens implant (images with intraocular lens implants were processed in the same way as the rest of the images). The average age of the subjects (at the time of the procedure) is 42 years and 57% of such subjects are female. We retrieved all the available images of the selected patients that corresponded to the mentioned scans. We had to discard around 1400 images because they didn't depict the whole anterior segment (they were deliberately zoomed in) and around 50 due to very poor quality. From those landmarks the following features were calculated: • Anterior chamber depth (L).
• Angle recess area (ARA_500, see 21 ) It has to be noted that some landmarks were not always clearly visible, in such cases the human expert guessed its position based on nearby features. The most common landmark that had to be guessed was the point in the top edge of the lens, which can be guessed based on the position of the iris. The sclerar spur is also frequently difficult to find but usually the clues were enough to guess its position by the expert. Ultimately, if the expert thought that it wasn't possible to make a good guess, the image was simply skipped and omitted from the subset of manually annotated images. As with the classification, these landmarks were not used by the algorithm, rather they were used afterwards to evaluate its results.  Unsupervised ordering and classification algorithm. In this section we present the algorithm for the unsupervised ordering and classification of OCT images. The input to the algorithm is a database of anterior chamber OCT images and the output is a map in a two dimensional plane. The algorithm performs three main steps: pre-processes the images, calculates a pair-wise distance measure between images and applies nonlinear dimensionality reduction.
Image pre-processing. The pre-processing of the images consists of the following substeps: homogenization, filtering, centering and aligning.  Homogenization: in each image, the intensity of each pixel was converted to double precision and normalized to be in between 0 and 1 (by linearly rescaling). Then, the horizontal and vertical spatial resolutions were adjusted to be the same (note that the original spacial resolution is anisotropic).   Filtering: First, a two dimensional rectangular median filter was applied to the image (with a 0.055 mm-by-0.117 mm rectangle). This filter was needed because of the process of adjusting the spatial resolution, which resulted in the noise being spread more in one direction than in the other. Then, an anisotropic diffusion 25,26 filter was applied to smooth the image, removing the spatial high-frequencies while preserving relevant edges. An example of such filtering is shown in Fig. 9.
For centering and aligning a set of statistical properties of the images were calculated, namely: where M(i, j) is the (gray) value of the image on the pixel that is i pixels down from the top edge and j pixels right from the left edge. With these quantities the centroid of the image (whose coordinates are ) and the covariance matrix (COV) were calculated:  Pair-wise distance measure between images. In this step a distance matrix (D(l, m)) was calculated whose entries are the pair-wise distances between images l and m. Two distance measures were employed: the Hellinger and Euclidean distances (d H ) 27 , Nonlinear dimensionality reduction. In order to extract meaningful information from the pair-wise distance matrix, D(l, m), we applied a nonlinear dimensionality reduction algorithm, directly to the pair-wise distance matrix. Two algorithms were tested: IsoMap 28 and t-SNE 29 .  These algorithms assign, to each image in the database, a point in a n-dimensional space, whose coordinates will be referred to as mapped coordinates. While the algorithm also works with an arbitrary number of dimensions, in this paper a two dimensional space (w, v) is used for visualization reasons. The choice of a two dimensional space is appropriated because the residual variance (defined as in 28 ) using two dimensions is of about 30%.
Computational runtime. All the described algorithms were implemented and run using MatLab in a portable computer with an Intel i7-7700HQ processor and 16 GB of RAM. We used the implementation of the nonlinear dimensionality reduction techniques written by van der Maaten et al. 30 . With this setup, it takes 5294 seconds (one and a half hours) to preprocess all the (1213) images including aligning and filtering, 1054 seconds (18 minutes) to compute the Hellinger distance matrix (735078 pair-wise distances), it takes 41 seconds for IsoMap to compute the mapping, and it takes 25 seconds for t-SNE to compute the mapping. It must be noted that all this runtimes could be significantly improved by rewriting the algorithms in a compiled language.
Informed consent statement. Informed consent statements were obtained from all the participants of the study.
Guidelines and regulations. OCT image acquisition and analysis were performed in accordance with the relevant European guidelines and regulations.