Three-Dimensional Spectral-Domain Optical Coherence Tomography Data Analysis for Glaucoma Detection

Purpose To develop a new three-dimensional (3D) spectral-domain optical coherence tomography (SD-OCT) data analysis method using a machine learning technique based on variable-size super pixel segmentation that efficiently utilizes full 3D dataset to improve the discrimination between early glaucomatous and healthy eyes. Methods 192 eyes of 96 subjects (44 healthy, 59 glaucoma suspect and 89 glaucomatous eyes) were scanned with SD-OCT. Each SD-OCT cube dataset was first converted into 2D feature map based on retinal nerve fiber layer (RNFL) segmentation and then divided into various number of super pixels. Unlike the conventional super pixel having a fixed number of points, this newly developed variable-size super pixel is defined as a cluster of homogeneous adjacent pixels with variable size, shape and number. Features of super pixel map were extracted and used as inputs to machine classifier (LogitBoost adaptive boosting) to automatically identify diseased eyes. For discriminating performance assessment, area under the curve (AUC) of the receiver operating characteristics of the machine classifier outputs were compared with the conventional circumpapillary RNFL (cpRNFL) thickness measurements. Results The super pixel analysis showed statistically significantly higher AUC than the cpRNFL (0.855 vs. 0.707, respectively, p = 0.031, Jackknife test) when glaucoma suspects were discriminated from healthy, while no significant difference was found when confirmed glaucoma eyes were discriminated from healthy eyes. Conclusions A novel 3D OCT analysis technique performed at least as well as the cpRNFL in glaucoma discrimination and even better at glaucoma suspect discrimination. This new method has the potential to improve early detection of glaucomatous damage.


Introduction
Optical coherence tomography (OCT) is a rapidly evolving technology and has gained a significant clinical impact in ophthalmology. [1][2][3][4] One of the major OCT applications is glaucoma assessment by measuring retinal nerve fiber layer (RNFL) thickness in circumpapillary and macula regions. Many studies showed that the current OCT quantitative assessment has excellent glaucoma discriminating ability. [1,3,5,6].
Spectral-domain (SD-) OCT's fast scanning speed allows threedimensional (3D) volume scanning of the retina, which may offer detailed and accurate quantitative analysis of the retinal structure. However, despite the rich information embedded in 3D OCT images, current standard quantitative structural OCT measurement is mostly limited to several hundred sampling points along a 3.4 mm circle diameter centered at optic nerve head, which does not take full advantage of the 3D dataset (over 20,000 sampling points). This sampling pattern was chosen mostly to allow compatibility with legacy data obtained using time-domain (TD-) OCT. The limited tissue sampling might lead to situations where early signs of structural changes are not detected when located outside the sampled circle (Fig. 1).
In addition to the circumpapillary RNFL (cpRNFL) thickness measurements, most SD-OCT devices provide the color coded thickness deviation map to highlight locations of structural damage (Fig. 1). These maps are generated by avaraging the mean RNFL thickness within a fixed number of neighboring sampling points (fixed-number super pixels) and comparing this measurement with age-matched normative databse. However, there is no quantitative summary of this super pixel analysis, and therefore clinicians' subjective interpretation is required.
In the clinical practice, subjects are often classified into three major groups: healthy, glaucoma suspect or glaucoma. Classifying subjects as ''suspect'' has important clinical advantage because some of these subjects posses pre-perimetric glaucoma characteristics. With the paucity of functional indication of glaucoma, correctly identifying glaucoma suspect based on their structural features is of upmost importance as it would allow timely adjustemt of clinical management. However, detection of preperimetric disease is still posing a significant challenge. [7,8] We hypothesize that a comprehensive use of the full 3D OCT data would improve detection of early glaucoma. The purpose of this study was to develop a novel 3D SD-OCT data analysis technique utilizing the full 3D dataset to improve the ability to detect glaucomatous structural damage at early stages. The new method uses variable-size super pixel mapping with a machine classifier analysis to quantitatively assess the full 3D OCT data. Unlike the conventional super-pixel analysis, which uses a fixed number of points, the newly developed variable-size super pixel is defined as a cluster of homogeneous adjacent pixels with variable size, shape and number. In this paper, we investigate whether the variablesize super pixel analysis is better suited for handling individual eye characteristics that might lead to improved glaucoma diagnosis predominantely in early disease stage.

Ethics Statement
This study was obtained the Institutional Review Board (IRB) approval named ''Optical Coherence Domain Reflectometry and Optical Coherence Tomography Measurements of Intraocular Structure''.

Subjects and Image Acquisition
This study was conducted with healthy, glaucoma suspect and glaucomatous eyes with a wide range of disease severity selected from the Pittsburgh Imaging Technology Trial (PITT). One hundred and ninety-two eyes of 96 subjects (44 healthy, 59 glaucoma suspect and 89 glaucomatous eyes) were enrolled to test the glaucoma discrimination performance. An independent dataset, including 46 eyes of 46 healthy subjects (randomly selected one eye for each subject), was used as the normative database. Institutional Review Board (IRB) approval was obtained for the study and all participants gave their consent to participate in the study. The study adhered to the Declaration of Helsinki and Health Insurance Portability and Accountability Act regulations.
Exclusion criteria for the study included history of ocular trauma or surgery (except for uncomplicated cataract or glaucoma surgery) and disease or treatment that might affect the visual field (e.g., stroke, central nervous sytem tumors) or retinal thickness (diabetes melitus, chronic steroid treatment).
All subjects had a comprehensive ophthalmic evaluation, reliable visual field (VF) and 3D SD-OCT scan all acquired within 6 months of each other. The ophthalmic evaluation included medical history, best-corrected visual acuity, manifest refraction, intraocular pressure (IOP) measurement, gonioscopy, slit-lamp examination before and after pupil dilation, and VF testing. VF was considered reliable if false negatives, false positives and fixation losses were less than 30%. All subjects had 3D OCT scans centered at the optic nerve head (ONH) (Cirrus HD-OCT; Carl Zeiss Meditec, Inc., Dublin, CA; ONH cube 2006200 scan protocol).

Clinical Diagnosis
Eyes were defined as healthy if there was no history of glaucoma, IOP#21 mmHg, the ONH did not meet the criteria for glaucomatous optic neuropathy (GON), as described below, and VF was normal.
Eyes were defined as glaucomatous if there was both GON and glaucomatous VF loss. GON was defined if either one of the following criteria was met: inter eye cup-to-disc (C/D) ratio asymmetry .0.2, accounting for disc size; rim thinning or notching; cup to disc ratio $0.6; RNFL wedge defect or disc hemorrhage. Glaucomatous VF was diagnosed if any of the following findings were evident on two consecutive VF tests: a glaucoma hemifield test outside normal limits, pattern standard deviation (PSD) ,5%, or a cluster of three or more non-edge points in typical glaucomatous locations (arcuate scotoma, nasal step, paracentral scotoma or temporal wedge), all depressed on the pattern deviation plot at a level of p,0.05, with at least one point in the cluster depressed at a level of p,0.01.
Glaucoma suspect eyes had IOP between 22 to 30 mmHg and/ or GON features all in the presence of a normal VF test.

3D SD-OCT Data Analysis
Super pixel mapping is a partitioning of image into a number of close-to-homogenous segments with variable sizes and shapes. The features of these variable-size super pixels were then extracted and used as the input for machine classifier analysis in order to generate a key index output for each image. The procedure to generate the index output included the following steps: comparison with normative database, feature map generation, super pixel mapping, feature extraction, and classification.
1. Normative Database. A normative database was assembled to measure the RNFL thickness deviation at each sampling point in the 3D dataset. Retinal layer segmentation, which was optimized for 3D SD-OCT dataset [9] was applied on each 3D OCT image to obtain the RNFL thickness and reflectivity at every single sampling point (total 2006200 points). All segmented 3D OCT datasets were visually evaluated to ensure the correct segmentation. Any OCT scan with .8% consecutive B-scan from the total number of B-scans that showed segmentaion error or .12% cumulative segmentation error was excluded from the study. ONH margin was automatically detected on the OCT fundus image using a software program of our own design. [10] Major retinal nerve fiber bundle path at each hemi-field on each RNFL thickness map was automatically detected. Each RNFL thickness map was then normalized by aligning the bundle location at each hemi-field to a reference position (population's average bundle location) for the comparison with the normative database in order to minimize spatial variability in RNFL thickness, especially at superior temporal and inferior temporal regions (Fig. 2). The RNFL thickness map normalization was processed on the concentrated circles with different radii started from the ONH center. For each subject, the ONH center was first aligned at a reference center point. At each concentrated circle, two average bundle location (superior and inferior) were computed from the normative database. Each subject's RNFL thickness profile at the given circle was normalized by strenching/shrinking so the subject's bundle location would coincide with the population average bundle location. The entire RNFL thickness      with the blood vessel mapping were extracted from each image (Fig. 3). The normalized RNFL map after applying the bundle path correction were compared with the normative database point by point to obtain a deviation map, with the cut-off value set to mean value minus stardard deviation (SD). Unlike the conventional setting of cut-off value, i.e., mean value -2SD in most OCT devices, this new setting is more sensitive to identify the case at the border of the normative database, which might be an indicator of the structural changes at the early stage. With this new setting, the RNFL data lower than the bottom 15.9% of the normative database was set with higher probability for the further process, comparing with the bottom 2.3% using the conventional setting. The internal reflectivity of retinal nerve fiber layer has been proved to be useful in glaucoma assessment. [11] The RNFL internal reflectivity, was calculated by taking the average reflectivity within the RNFL along each A-scan, with each voxel's reflectivity normalized to its A-scan's saturation before the averaging. Retinal blood vessel generated shadow at retina nerve fiber layer, which was noise in the RNFL thickness computation. Moreover, the vessel patterns varied randomly among subjects. To minimize the vessel effect on the RNFL thickness map, the retinal blood vessels in the 3D dataset were automatically detected using a 3D boosting algorithm [12] and filled out. The RNFL thickness of each pixel located at the blood vessel region was replaced by a value computed from all the non-vessel pixels on the RNFL thickness map using bi-linear interpolation. The final feature map is a function of the RNFL thickness and interal refleciviy after accounting for the blood vessls and the deviation map.
3. Variable-size super pixel segmentation. Variable size/ shape super pixels were automatically mapped on the 2D feature map by grouping homogeneous neighboring pixels using a ncut algorithm. [13] The ncut algorithm is to partition an image into dozens to thousands of small regions (called super pixels) by grouping neighboring pixels, where pixels within a partitioned region have homogenous properties while different partitioned regions have maximal differences in their properties. One hundred super pixels were initially segmented on the feature map. The size of each super pixel was automatically adjusted with the predefined criteria based on the pathologic contexts of glaucoma. To be more sensitive to RNFL thinning (glaucomatous damage), smaller super pixels were assigned to thinner RNFL. Each initially segmented super pixel was recursively partitioned into N more super pixels, while N was a function of mean, standard deviation, size, and deviation to the normative database of the given super pixel compared with global mean and standard deviation of the 2D feature map, written as: if h(Z sp ,s sp ,N sp ,Dev sp )~0 and Z sp vZ avg where Z sp ,s sp ,N sp ,Dev sp represented the average RNFL thickness, RNFL standard deviation, size and average deviation of the given super pixel, N min was a pre-defined minimal size of super pixel, Z avg and Z max were the average and maximal thickness with the entire 2D RNFL thickness map. h was defined with a series of criteria with ''IF'', ''logical AND'', ''logical OR'' and ''logical NOT'' operations, which corresponded to different super pixel conditions need further partition or not, written as: a 1 and a 2 were two constants chosen based on the segmentation logic, which were used to control the number of segments partitioned in the given super pixel. For a 1 , because super pixel RNFL thickness was thicker than the average RNFL thickness, this super pixel was further partitioned only if its standard deviation was large enough. We wanted the super pixel partitions into several big segments, and recursively partitions into small segments in the following iterations. Therefore, a 1 was set to 0.2 to lower the number of the segments. For a 2 , since the super pixel RNFL thickness was thinner than the average RNFL thickness, we wanted to partition the super pixel into many small segments. Therefore, a 2 was set to 1.2 to higher the number of segments. The stability is an important parameter in ncut algorithm to control the super pixel segmentation result. It was set to a small value, i.e., 20.1, to let the 100 initial super pixels have more flexible boundaries. In the recursive partition processing, to make In summary, the super pixel number and size were automatically adjusted by initial segmentation and recursive partition. The segmented feature map provided a qualitative analysis with more natural representation, in which damaged areas tended to have smaller super pixels, while normal regions had larger super pixels (Fig. 4).
4. Feature Extraction. Super pixel map provided a qualitative representation of the structural damage. To summarize the map as quantitative disease indices, a total of 68 super pixel features were extracted and used as the inputs of machine learning classifiers analysis (Table 1). Three thresholds, Tmin, T1, and T2 were set to 50, 133, and 400 pixels respecitvely based on the experiments to obtain two sub-groups of super pixel with large size and small size. The features were calculated from all super pixels together with two sub-groups.
5. Glaucoma Classification. Glaucoma classification was performed by implementing LogitBoost adaptive boosting algorithm, [14] which was designed as a supervised two-class machine classifier. Because the machine classifier was a two-class classifier, only two of the three clinical groups (healthy, glaucoma suspects, and glaucoma) were used to train the classifier each time. The output of the classifier was a continuous number ranging from negative (health) to positive (disease) value. Ten-fold cross validation was used to train and test the machine classifier.

Statistical Analysis and Evaluation
The software performance was evaluated by the area under the receiver operating characteristics (AUCs) of the machine classifier outputs, for discriminating between healthy and diseased eyes.
AUCs were compared with those of the conventional method of diagnosis -the mean and best quadrant measurements of the circumpapillary RNFL (cpRNFL) thickness generated by Cirrus HD-OCT software. If both eyes had the same clinical diagnosis, one eye was randomly selected from each subject to compute AUCs. AUCs were compared using the Jackknife method, [15] which is equivalent to the method of DeLong et al. [16] Sensitivity at the 85% specificity was calculated for each of the machine classifier outputs and compared with the measurements of conventional cpRNFL analysis.

Results
One hundred and ninety-two eyes of 96 subjects (44 healthy, 59 glaucoma suspect and 89 glaucomatous eyes) qualified for the study. The characteristics of the study participants were summarized in Table 2. Five glaucomatous eyes (5 subjects; 5.2%) were excluded from the study due to retinal layer segmentation failure.
Examples of the results of super pixel mapping are given in Fig. 4. The super pixel boundaries were superimposed on the processed grey scale RNFL thickness map, where a brighter pixel intensity corresponded to a thicker RNFL. Different super pixel distribution patterns are evident in the various stages of glaucomatous damage. The population's average distributions of super pixel RNFL thickness and size for healthy, glaucoma suspects, and glaucomatous eyes are illustrated in Fig. 5.
Glaucoma discriminating performance was assessed with three different grouping combinations: healthy vs glaucoma+glaucoma suspects (HvGGS), healthy vs glaucoma suspects (HvGS), and healthy vs glaucoma (HvG). Table 3 gives the performance of cpRNFL thickness measurements in global and 4 quadrants. Comparing with the conventional cpRNFL thickness, machine classifier provided better AUCs and higher sensitivities (Table 4 and Fig. 6). The AUC for HvGS was statistically significantly higher with the super pixels analysis than the conventonal cpRNFL thickness (0.855 vs 0.707, respectively, p = 0.031, Jackknife test), while no significant difference was detected with HvGGS and HvG. Comparing with the best quadrant measurements, i.e., the inferior quadrant, the machine classifier showed higher sensitivities for HvGGS and HvGS without reaching statistical significance (Table 5).

Discussion
Current standard quantitative glaucoma analysis using 3D SD-OCT have substantial room for improvement especially for detection of early glaucoma. [7,8] One possible reason is that only a small fraction of the 3D dataset is used in the analysis. We suggest a new 3D OCT data analysis method, using machine  learning technique based on variable-size super pixel mapping to quantitatively summarize the full 3D dataset and automatically identify glaucomatous eyes. We demonstrated that this method was better at discriminating glaucoma suspect eyes from healthy eyes (HvGS) comparing with mean cpRNFL thickness, and performed at least as well for HvG and HvGGS. The improved performance in discriminating early disease is because this method is tailored to detect localized damage that typically occurs at early stage of glaucoma. However, in late stage of the disease, a globalize damage is more common where the variably sized super pixel analysis has the same performance as the cpRNFL analysis.
There are several advantages of variable size/shape super pixel machine classifier analysis. First, flexible size and shape super pixel processing provides more natural representation to fit the variable spatial architecture of the structural damages. Many ocular diseases demonstrate areas of pathologic change with measureable differences from unaffected areas when considering features such as retinal layer thickness, internal reflectivity, etc. [6,[17][18][19] These pathologically affected areas usually share similar characteristics but with variable magnitude in variable shape and size. Therefore, variable-size super pixel segmentation efficiently depict this natural representation by grouping similar neighboring sampled points based on the homogeneity of various features. Second, super pixel processing is computationally efficient. Super pixel mapping reduces the number of data points from over 20,000 sampling points down to a few hundred super pixels by condensing the redundant information. This is a common approach for handling such a high density data with the balance of preserving meaningful information and improving the computational efficiency. Moreover, machine classifier provides an efficient and optimized way to combine dozens of super pixel features together with global RNFL features into one key index to automatically identify diseased eyes. The performance of the software was based on the combination of super pixel processing and machine learning technique. Our study showed that the new super pixel analysis method quantitatively summarize the full 3D OCT dataset, which outperformed the conventional cpRNFL analysis in terms of discriminating ability and diagnostic sensitivity for early glaucoma detection.
To make this super pixel analysis sensitive to the localized damages, we set a criterion of super pixel mapping so that smaller super pixel corresponded to the thinner RNFL. With the natural distribution of the retinal nerve fiber, the RNFL is thicker in superior-temporal and inferior-temporal areas and thinner in temporal and nasal areas. In our method, the deviation from the normative database was used to partially control the super pixel mapping (Equation 1 and 2). Therefore, although temporal and nasal areas had relatively thinner RNFL thickness, they would not have smaller super pixels if the RNFL thickness was within the range of the normative database.
In the glaucoma classification step, both local features (super pixel feautres) and glocal features (global RNFL measurements) were used to train and test the boosting machine classifier. Global featuares, such as cpRNFL, were able to discriminate most glaucomatous eyes. Local features were able to enhance localized glaucomatous damages for eyes at the early stages of the disease. Therefore, combining both glocal and local featuares in the machine classifier improved the diagnostic sensitivity comparing with the conventional cpRNFL measurement. Including more features and applying feature selection operation may further enhance the performance of the machine classifier.
Searching the literature, we were able to find only one publication where glaucoma strctural analysis used the information from the entire 3D dataset. [7] In that study, 2D RNFL deviation map was analyzed and subjectively graded into 5 different groups. Compared with conventional cpRNFL analysis, this analysis provided additional spatial and morphologic information of RNFL damage, and significantly improved the diagnostic sensitivity for glaucoma detection. This is consistent with the result of our super pixel analysis study that analyzed the full 3D dataset providing more localized and detailed information of structural damages, and showing the potential to detect structural damage in early stages of the disease. The advantages of our method are the fully automated process and advanced data analysis compared to the subjective grading system used in the previous study. The expansion of this analysis will be including the distribution, shape and spatial locations of the localized structural changes with the super pixel and machine classifier analysis, which may further improve the performance of the 3D data analysis.
In this study we used only one SD-OCT device out of a variety of devices that are commercially available. However, the conceptual approach can be applied to any of these devices.
In conclusion, the super pixel processing with machine classifier analysis generated from 3D SD-OCT data has the potential to improve early detection of glaucomatous structural damages. This method can be easily extended to other ocular diseases by modifying the features corresponding to the various pathologic contexts.