Analysis of GDx-VCC Polarimetry Data by Wavelet-Fourier Analysis across Glaucoma Stages

Abstract

purpose. The purpose of this study was to apply shape-based analysis techniques of retinal nerve fiber layer (RNFL) thickness to GDx-VCC (variable corneal and lens compensator; Laser Diagnostic Technologies, Inc., San Diego, CA) polarimetry data and to evaluate the techniques’ ability to detect glaucoma in its earliest stages. Wavelet-based (wavelet-Fourier analysis [WFA]), Fourier-based (fast Fourier analysis [FFA]), and several previous variations of shape-based analysis were considered, as well as the standard metric nerve fiber indicator (NFI), and all were compared as a function of disease stage.

methods. GDx-VCC scans of one eye of each of 67 patients with glaucoma and each of 67 healthy age-matched subjects provided RNFL thickness estimates at a fixed distance from the optic disc. Severity of disease was graded according to the Glaucoma Staging System and also by mean deviation (MD) from standard automated perimetry. WFA, FFA, and NFI procedures were performed including the following variations: use of signed or unsigned phase, inclusion of interocular or intraocular asymmetry of analysis parameters, and combination of features by principle components analysis or Wilks λ. Independent samples (k-fold variation) were used for training and testing. Sensitivity, specificity, and receiver operating characteristic (ROC) area were obtained.

results. Classification performance of WFA (ROC = 0.978) was significantly better than FFA (ROC = 0.938) and NFI (ROC = 0.900). This difference was largest for the earliest stages of glaucoma. Shape-based analysis methods performed better than NFI overall. Adding between-eye asymmetry measures helped FFA but not WFA.

conclusions. Shape-based analysis, and WFA in particular, makes an important improvement in detecting earliest glaucoma with polarimetry.

Glaucoma is a progressive optic nerve disease characterized by nerve head changes and characteristic visual field loss on standard perimetry. However, it has become apparent that this visual field loss in glaucoma occurs only after significant neural damage has occurred. ^{1 2 3 4 5 6 7 8 9} This conclusion makes it even more important to establish new tests of function and/or neuroanatomic loss that better delineate and characterize early glaucoma. To that end, the present paper reports our efforts to improve the characterization of neuroanatomic disruption in glaucoma. Specifically, we assessed the utility of inferring local thickness of the retinal nerve fiber layer (RNFL) by polarimetry with the GDx-VCC (variable corneal and lens compensator; Laser Diagnostic Technologies, Inc., San Diego, CA). 

At present, it is not known how best to organize and extract clinically useful information from the many (approximately 65,000) local thickness values provided by clinical polarization (or tomography) devices. One approach has been to characterize mathematically the shape of the surface conveyed by the map of these local thickness measurements (Essock EA, et al. IOVS 1999;40:ARVO Abstract 3481; Essock EA, et al. IOVS 2001;42:ARVO Abstract 93; Sinai MJ, et al., IOVS 2001;42:ARVO Abstract 717; Sinai MJ, et al. IOVS 2002;43:ARVO E-Abstract 302; Bryant FD, et al. IOVS 2000;41:ARVO Abstract 485; Essock EA, et al. IOVS 2003;44:ARVO E-Abstract 3378; Gunvant P, et al. IOVS 2004;45:ARVO E-Abstract 5504). ^{10 11} That is, if the shape of the data surface that is characteristic of healthy eyes can be quantified, then deviations from the normative shape could be quantified. In addition, emphasizing shape rather than thickness can obviate the large individual variation among healthy eyes of general thickness. ^{12 13 14} In the typical analysis, RNFL thickness measurements at points at a given distance from the disc as a function of angle are analyzed at one or more distances from the disc (i.e., one or more “rings”). This TSNIT (temporal, superior, nasal, inferior, temporal) graph of thickness for a ring conveys a general double-hump pattern of thickness due to the much greater number of ganglion cell axons entering the disc superiorly and inferiorly. However, the newer imaging techniques indicate that the smooth, “textbook” double-hump pattern of TSNIT RNFL thickness of an individual healthy eye can actually be quite jagged, often with very low local minima at various locations within the superior and inferior humps of the TSNIT pattern. In addition, it has been reported that a significant number of healthy eyes show a bimodal or “split bundle” (Vermeer KA, et al. IOVS 2004;45:ARVO E-Abstract 3309; Varma R, et al. IOVS 2002;43:ARVO E-Abstract 262). ¹⁵  

Numerous mathematical methods could be used to analyze and parameterize the TSNIT shape. Previously, we used Fourier analysis as the mathematical analysis method to characterize the shape of the TSNIT pattern (Fig. 1) . This “fast” Fourier analysis (FFA) linearly breaks up the pattern into sinusoidal variations in thickness (i.e., into a set of sinusoids in which each sinusoid is a different scale, or frequency) and thus has a different number of humps across the TSNIT data set. ¹⁰ Fourier analysis breaks an individual’s TSNIT pattern into the sum of a set of sine-wave patterns of particular amplitude and phase (position). We have applied this analysis method to the TSNIT pattern previously in three ways: simultaneously at multiple distances from the disc in a two-dimensional polar Fourier analysis (Essock EA, et al. IOVS 1999;40:ARVO Abstract 3481) at a single ring with the Fourier analysis performed separately on the superior and inferior halves of the data (Essock EA, et al. IOVS 1999;40:ARVO Abstract 3481; Essock EA, et al. IOVS 2001;42:ARVO Abstract 93; Bryant FD, et al. IOVS 2000;41:ARVO Abstract 485; see also Ref. ¹⁶ ); and most typically, at a single ring for the full TSNIT pattern (Essock EA, et al. IOVS 2001;42:ARVO Abstract 93; Sinai MJ, et al. IOVS 2001;42:ARVO Abstract 717; Sinai MJ, et al. IOVS 2002;43:ARVO E-Abstract 302; see also Ref. ¹¹ ). ¹⁷ We have also considered combinations of amplitudes of specific frequencies (i.e., specific shapes) (Essock EA, et al. IOVS 1999;40:ARVO Abstract 3481; Essock EA, et al. IOVS 2001;42:ARVO Abstract 93; Bryant FD, et al. IOVS 2000;41:ARVO Abstract 485) ¹⁰ both unsigned phase (i.e., indicating magnitude) (Sinai MJ, et al. IOVS 2001;42:ARVO Abstract 717; Sinai MJ, et al. IOVS 2002;43:ARVO E-Abstract 302) ¹⁷ and signed phase (indicating magnitude and direction) (Essock EA, et al. IOVS 2001;42:ARVO Abstract 93) and measures of asymmetry of these amplitude and phase parameters, both between the fellow eyes and also between superior and inferior hemiretinas within an eye (Essock EA, et al. IOVS 1999;40:ARVO Abstract 3481; Essock EA, et al. IOVS 2001;42:ARVO Abstract 93; Sinai MJ, et al. IOVS 2002;43:ARVO E-Abstract 302; Bryant FD, et al. IOVS 2000;41:ARVO Abstract 485). The use of asymmetry values has appeared to be a particularly promising addition to FFA and are directly examined herein. However, despite the success of these various approaches, a Fourier analysis may not be the best choice to capture the shape of the normal TSNIT pattern in relatively few parameters. In the present paper, we report use of a different shape analysis method (Essock EA, et al. IOVS 2003;44:ARVO E-Abstract 3378; Gunvant P, et al. IOVS 2004;45:ARVO E-Abstract 5504) ^{18 19 20} that was developed to address this weakness of using FFA for RNFL analysis. Instead of a Fourier analysis that emphasizes the frequency domain (i.e., a set of the infinite regular repeating patterns of sine waves), we adopted a wavelet analysis for the primary analysis of the TSNIT pattern. A wavelet analysis offers the advantage that it emphasizes local shape (i.e., the space domain) and is better suited to capture the irregular or abrupt changes in the TSNIT shape that are particularly evident in GDx-VCC data (see Fig. 6B ) (Vermeer KA, et al. IOVS 2004;45:ARVO E-Abstract 3309). ²¹  

The first objective in the present study was to compare this new wavelet-based TSNIT shape analysis method, wavelet-Fourier analysis (WFA), to several variations of the previous Fourier-based analysis method and to the standard GDx-VCC metric provided by the device (nerve fiber indicator, NFI). The second objective was to consider the ability of these TSNIT data analysis methods to detect glaucoma as a function of disease stage; as, to be useful clinically, RNFL analysis must be able to discriminate earliest glaucoma from the RNFL pattern in healthy eyes. 

Methods

Subjects

One-hundred and thirty-four eyes (70 OD and 64 OS) of 134 individuals (67 healthy subjects and 67 with glaucoma) were included prospectively from the outpatient clinics of glaucoma specialists at The Institute of Ophthalmology and Visual Science, University of Medicine and Dentistry of New Jersey (UMDNJ) (New Jersey Medical School, Newark, NJ); the New York Eye and Ear Infirmary (New York, NY); and the Eye Care Center (San Diego, CA), and the data were provided by Laser Diagnostic Technologies, Inc. for this purpose. Data management and data analysis conformed to HIPAA (Health Insurance Portability and Accountability Act) regulations, institutional review board approval was obtained, and the protocol adhered to the provisions of the Declaration of Helsinki. The groups were matched closely for age with the mean age of the healthy-subjects group being 64.61 and 67.22 years for the glaucoma group. The difference in age between the groups was not significant (P > 0.05). None of the study participants had significant ocular media opacity that could affect their visual field results or imaging with GDx-VCC. 

Visual Fields

Automated static threshold perimetry was performed (Humphrey Field Analyzer II; Carl Zeiss Meditec, Inc., Dublin, CA). All visual fields had good reliability and were either a central 30-2 or a 24-2 threshold test with size III white stimulus with one of two algorithms (standard or full Swedish Interactive Threshold Algorithm). The severity of visual field defects was analyzed by staging the visual fields on the basis of the Glaucoma Staging System [GSS] ^{22 23} and were also staged with respect to their mean deviation (MD). All healthy subjects had normal intraocular pressure and had normal appearance of optic discs. Visual fields were measured in most, but not all, cases and were normal. 

Full details of the GSS are provided elsewhere. ²² Briefly, it is a nomogram that uses the MD and pattern standard deviation (PSD) on a Cartesian coordinate diagram with MD on the x-axis and corrected PSD on the y-axis. The GSS can be used to stage visual field damage in any field report that gives the visual field indices (MD and the corrected PSD). If the corrected PSD was not available or was unreliable, the corrected PSD was calculated by adding 0.7 to the PSD. ²² The GSS classifies subjects into six categories: normal, and stages 1 to 5 of severity. ²² All patients in the glaucoma group had a visual field defect of GSS stage 1 or greater. The stages 1 to 3 are considered to be mild to moderate severity with stage 1 considered to have very subtle defects; stage 2 is usually field plots with nasal steps or mild scotoma; and stage 3 is usually an overt defect such as dense arcuate scotoma. The visual field at stage 4 is considered to be at an advanced stage of disease and a stage 5 visual field has a very low threshold with only small remnants of sensitivity. The GSS can also classify the visual field defects into three other categories: generalized defects, localized defects, and mixed defects. Generally, the mixed defect is observed to be the most common type of defect, whereas patients with early glaucoma are typically observed to have localized defects. ²² The present sample was consistent with this finding. 

Average MD of the visual field of all patients was −6.82 ± 6.2 dB (SD), and average PSD was 6.25 ± 4.2 dB (SD). The patients were grouped into five stages based on their MD: stage 1 (0 to −5 MD), stage 2 (−5 to −10 MD), stage 3 (−10 to −15 MD), stage 4 (−15 to −20 MD), and stage 5 (−20 to −25 MD). The stage of the disease according to GSS and the MD were used to identify patients with early field defect and to examine further the performance of the WFA in different stages of disease. 

Scanning Laser Polarimetry

The GDx-VCC was used for this study. The measurements were obtained in three different clinics and were performed by experienced technicians. The GDx-VCC calculates and compensates for the anterior segment birefringence for each eye using the method described by Zhou and Weinreb. ²⁴ Macular images are obtained, and the combined polarization of the radially oriented axons of Henle’s fiber layer of the macula and the anterior segment (a bow-tie pattern on the retardation image) is analyzed. A software program within the GDx-VCC then determines an estimate of polarization magnitude and axis that is due to the anterior segment birefringence. These estimates are applied to compensate a polarization retardation map at the disc for the individual eye’s anterior segment contribution to such maps. Retardation is converted to an estimate of depth ²⁵ and the software provides 256 RNFL depth estimates within an eight-pixel-wide ring with inner and outer radii of 27 and 35 pixels, respectively. These estimates are grouped into 64 sectors to yield a TSNIT curve of 64 points in the 360° around the disc. 

Application of Shape-Analysis Methods

Wavelet-Fourier Analysis.

There are three major steps involved in wavelet-Fourier analysis (WFA): feature extraction, feature optimization, and classification. First, features are extracted from the GDx data. We used a discrete wavelet transform (DWT), which is a common method used in feature extraction ²⁶ and is highly suitable for analyzing discontinuities and abrupt changes contained in signals. ²⁷ The wavelet transform applies the scaling function (Fig. 2A)and the wavelet function (Fig. 2B)yielding a down-sampled signal represented by two sets of coefficients, termed the approximation part and the detail part. The approximation coefficients contain down-sampled spatial information, and the detail coefficients contain detailed information characterizing the difference between the signal and the approximation (Fig. 3) . DWT is a multiscale analysis method, meaning that the analysis can be based on various scales. ²⁷ At each transformation scale (or “level” of this multiscale pyramid), there is a different resolution ability encoded in that level’s approximation and detail coefficients (i.e., at each level, the DWT analysis is applied to the prior level’s approximation coefficients). We found that a second-level DWT analysis works well for 64-point VCC data. The second-level decomposition preserves high enough spatial resolution ability (in the approximation part), but results in a relatively low-frequency resolution of the detail part. For this reason, a further step, drawing from our previous work, ¹⁰ applies a Fourier transform to the second-level DWT detail part (constituting a second-step transformation), and provides frequency amplitudes to achieve high frequency resolution for WFA. The final part of this step is to normalize DWT approximation coefficients and FFT amplitudes separately, and to join them together as the final set of features used in the analysis (i.e., the “feature vector”). Second, the dimensions of the feature vectors are reduced by using principal components analysis (PCA) ²⁷ to combine (and weight) various features to maximize the spread of the data points in the resultant multidimensional (reduced) feature space. The resultant lower dimensionality of the feature vectors can make the classifier more efficient and more stable. Third, a classification procedure is applied to the data. The Fisher linear discriminant function (LDF) was used as a classifier, as it performs slightly better than other similar methods. ²⁸ The role of an LDF is to provide a criterion that optimally classifies a set of values into two categories (in the present study: glaucoma and healthy). To assure external validity of the obtained LDF, the LDF must be tested on data (a “test” set) that are independent of the “training” set from which they were derived. To achieve this, the dataset was first randomly split into two independent subsets, one for training the classifier and obtaining the LDF, and one for testing it. To provide a robust validation procedure, a variation of k-fold cross validation ^{29 30} was used, as it is superior to methods such as split-half or holdout and is especially advantageous when a very large sample is not possible. The procedure is similar to the split-half procedure, but the estimates are based on the mean of estimates from k random splits, where k is usually several hundred. This method reduces bias and increases reliability. ^{29 30} Also, this repetition allows the use of a larger proportion of the cases for the training set, resulting in greater stability of the final estimate. The method has been used effectively previously in similar studies, ^{20 31} and, in effect, is a resampling procedure comparable to the bootstrap method. ³² For this procedure, the majority of cases (in our study, 90%) are selected for the training set, and the classification is applied to the small test set consisting of the remaining cases (10%). Then a new random selection is obtained and classification performed, and this procedure is repeated k times (in our study, 400 times) and the results averaged. 

For the particular wavelet pattern, we chose a fourth order of Symlets (Fig. 2) , although numerous wavelets are likely to work. Symlets is a modification of Daubechies wavelets, the symmetry of which is increased while retaining great simplicity. First, we applied the two-scale DWT to the 64-sector GDx-VCC data, then retained the second-level approximation coefficients but applied an FFT to the second-level DWT detail coefficients (Fig. 3) . Note that we used only the coefficients of the second scale and discarded those of the first scale. We found that the coefficients of the two scales were so similar that they would be suppressed by the PCA even if both were used as features. Coefficients of only the amplitudes, and no phases, from the FFT were used. Next, the DWT coefficients and FFT coefficients were each normalized to the range (0–1) and were joined as preliminary feature vectors. It is difficult to decide the number of the reduced features to use (n) in PCA. We varied n and empirically evaluated the classification performance (using area under the ROC curve to compare performance) which is shown in Figure 4 . The best value of n was found to be 4 (where a global maximum of the area under the curve was achieved; Fig. 4 ) and this value was selected for all subsequent PCA and classification. These optimized features with lower dimensionality were then supplied to the Fisher linear discriminant analysis (LDA) program, which produced an LDF based on the training set. Finally, this LDF was applied to classify the independent test set, and the parameters used to evaluate performance (sensitivity/specificity and area under the curve [AUC]) were obtained on 400 runs and averaged. All results for FFA and WFA reported herein, except as noted, are based on these averages of 400 runs. 

Fast Fourier Analysis.

The FFA, performed to compare with WFA performance, was conducted as it has been previously (Essock EA, et al. IOVS 1999;40:ARVO Abstract 3481; Essock EA, et al. IOVS 2001;42:ARVO Abstract 93; Sinai MJ, et al. IOVS 2001;42:ARVO Abstract 717; Sinai MJ, et al. IOVS 2002;43:ARVO E-Abstract 302; Bryant FD, et al. IOVS 2000;41:ARVO Abstract 485; Essock EA, et al. IOVS 2003;44:ARVO Abstract 3378). ^{10 11 17} The basic FFA method is to perform an FFT on the TSNIT data set (in this study, 64 points). All amplitude and phase coefficients (and the DC, or offset value) are entered into a procedure to obtain a classification function. In the present work, we compared the utility of several variations of this method, most of which have appeared previously in the literature. The following comparisons were made: The FFA features were optimized either by Wilks λ, which has been the main method used previously for FFA (Essock EA, et al. IOVS 2001;42:ARVO Abstract 93; Sinai MJ, et al. IOVS 2001;42:ARVO Abstract 717; Sinai MJ, et al. IOVS 2002;43:ARVO E-Abstract 302), ^{11 17} or by PCA, as used in the WFA procedure. The phase coefficients from the FFT were coded as either signed or unsigned (i.e., absolute values taken). Signed coefficients indicating the direction (clockwise or counterclockwise) as well as magnitude that a given sine-wave component was located on the TSNIT pattern compared with another component. Unsigned coefficients indicated only the magnitude of a shift of the given component relative to a reference. Asymmetry of the FFT values was added as another set of features, in addition to the amplitude and phase values produced by FFA. 

For between-eye asymmetry, an FFA was also completed on the fellow eye and the absolute arithmetic difference between a given FFA feature was obtained in the two eyes. Thus, for this analysis, the DC level of the test eye, the 32 amplitudes, the 32 phase values, and 65 asymmetry values (all amplitude, phase, and DC values) were used as the feature set (i.e., a total of 130 features). For WFA, the 32 features of the test eye and the 32 values obtained from the absolute arithmetic difference of these values for the fellow eyes were used as features. Hence, 64 features were used before application of PCA. For NFI, to add a measure of between-eye asymmetry, the absolute arithmetic difference between the NFI of the fellow eyes and the NFI of the test eye were combined in a Fisher LDA procedure to obtain the optimal classifier. 

For within-eye asymmetry measures, data from separate analyses of the superior 32 points and the inferior 32 points was compared by taking the absolute difference in transformed feature space. For FFA, the eye’s full 64-point transformed coefficients (65 features) were combined with 33 within-eye asymmetry features, creating a total of 98 features to be used as preliminary features fed to the PCA. Similarly, for WFA, 32 transformed features plus the 20 within-eye asymmetry features were used. For NFI, within-eye asymmetry was not possible to calculate without changing the nature of the NFI measure, so this measure was omitted. 

For comparison purposes throughout this article, FFA calculated with unsigned phases, without asymmetry measures, and using PCA was considered to be the “standard” FFA method. 

Statistical Analysis

Independent sample or correlated t-tests were used as appropriate with a significance level of 0.05 adopted. To evaluate the classification performance of the three metrics, WFA, FFA, and the NFI, the sensitivity, specificity, and the receiver operating characteristic (ROC) curves were calculated. ROC curves were generated in the standard fashion by systematically varying the cutoff value and noting performance. In the instances of numerous (400) runs (i.e., for the WFA and FFA methods) test cases from all runs were plotted to provide the most accurate estimate of the mean ROC curve. To compare the ROC curves, the following procedure was used. For each particular pair-wise comparison (e.g., WFA versus FFA), the difference in the areas of the ROC curves (based on the 400 k-fold samples) were then easily computed (ROC area for WFA minus ROC area for FFA). These could then be used to compute approximate confidence intervals, by using the percentile method. ³² Bonferroni-corrected confidence intervals, each with a confidence level of 98.3% were computed as 0.83% and 99.17% to achieve an overall coverage of 95%. ³³  

Results

Comparison of Performance of Analysis Methods

Both shape analysis methods were better able to classify glaucoma from healthy eyes than the standard GDx-VCC classifier NFI as shown in Table 1 . Furthermore, the new shape analysis method, WFA, considerably outperformed the FFA method with an ROC area of 0.978 compared with 0.938 for the FFA method (P < 0.0001), in addition to outperforming the NFI method (ROC area of 0.900; P < 0.003). (However, the obtained difference in FFA and NFI performance (0.938 vs. 0.900) was not statistically significant (P > 0.05), with FFA producing larger ROC areas in 92% of the samples.) A similar trend showing the best performance with WFA, worst with NFI, and intermediate with FFA, is also seen in the sensitivity at a fixed specificity (90% or 95%) as shown in Table 1 . For example, when specificity is fixed at 90% or 95%, sensitivity with WFA is 92.5%, and 87.8%, respectively, considerably higher than that with the other methods. 

Whether Wilks λ or PCA was used to evaluate the features of FFA made little difference overall (Tables 2 3) . Averaged across the six variations of FFA (that is, with either signed or unsigned phase, for each of the following three methods: with no asymmetry, with between-eye asymmetry, or with within-eye asymmetry), ROC AUC was quite similar: 0.929 for the Wilks λ version and 0.945 for the PCA version (P > 0.05, correlated t-test). There was no evidence of a systematic advantage of either PCA or Wilks λ (Table 2 3) . 

In general, whether adding asymmetry measures improved glaucoma classification performance depended on which analysis method and type of asymmetry measure were considered (Table 2) . Inclusion of between-eye asymmetry of Fourier phases and amplitudes improved classification performance. Considered across all Fourier methods (all combinations of PCA/Wilks and signed and unsigned phases), inclusion of between-eye asymmetry increased ROC area from 0.927 to 0.969 (P < 0.05, correlated t-test). Including asymmetry information in the NFI helped only slightly, but not significantly (ROC AUC of 0.900 increased to 0.911; P > 0.05). WFA performance, however, was not affected by inclusion of between-eye asymmetry of its features (ROC AUC, 0.978 and 0.977). Use of within-eye asymmetry did not help any measure and appeared to hurt performance at times. 

Whether phase was used as signed (i.e., identifying a clockwise or counterclockwise position shift) or unsigned (simply magnitude of offset) made very little difference in the classification results of the Fourier method, whether between- or within-eye asymmetry was used. Averaged over all methods, performance with unsigned phase was 0.938 and 0.937 for signed phase. Table 3shows these results and also the full results for all measures considered in this study. 

TSNIT Shape Characteristics Used to Discriminate Glaucomatous Eyes

The mean TSNIT thickness curve for the glaucomatous and healthy eyes is plotted in Figure 5D . The figure illustrates the shapes of the four features used in the WFA classification procedure. The characterization of healthy and glaucomatous eyes by only the first or only the second feature is shown in Figures 5A and 5B . The fit to the mean curves for the healthy and glaucoma groups if the first two features are used is shown in Figure 5C , and the final fit with all four features used is shown in Figure 5D . For comparison, Figure 5Dalso shows the means of the original data. Although all four features contribute to WFA performance, mean curves for either the third or fourth features showed very little difference in the glaucomatous and healthy eyes. Examples of the fit to individual eyes are shown in Figure 6 , for an eye with a smooth TSNIT pattern (Fig. 6A)and for an eye with a fairly jagged VCC TSNIT curve (Fig. 6B) . 

To compare the two shape analysis methods, Figure 7shows plots of the mean-thickness curve and both the mean-WFA curve and mean-FFA (PCA) curve fits to the data of all healthy (Fig. 7A)and glaucomatous (Fig. 7B)eyes. These curves differed in a variety of ways, but differences were most apparent at the peaks and troughs. Overall, the fit of the WFA curve is somewhat better than that of the FFA curve for the patients with glaucoma, and the fits are equivalent for the healthy eyes. The fits were quantified by calculating the shape deviation (σ_s) between the fitted and actual data:  

\[{\sigma}_{\mathrm{s}}\ {=}\ \sqrt{{{\sum}_{i{=}1}^{64}}(t{^\prime}_{i}\ {-}\ t_{i})^{2}/(64\ {-}\ 1)}\]

, where t′ and t _i are RNFL thickness for reconstructed and original curves, respectively, yielding a shape deviation of 3.5 and 5.3 for the WFA and FFA fits, respectively, in glaucomatous eyes and 4.8 and 5.0 in the healthy eyes. 

Classification Performance Stratified by Visual Field Defect

To truly have clinical utility, a classification method must do more than classify patients whose disease is otherwise apparent. To assess the ability of shape analysis to detect earliest glaucoma, we stratified the glaucomatous eyes by severity of visual field loss and compared the ability of the WFA, FFA, and NFI methods to correctly classify glaucomatous from healthy eyes in each category (Table 4) . GSS staging ²² is shown in the top part of the table and severity of MD loss is shown in the middle part. From Table 4 , it can be seen that it is the earlier/milder visual field losses in which WFA and NFI performances differ most. At the mildest levels of disease (i.e., stage 1 GSS and MD < −5.0), the performance of both shape-based methods is considerably superior to that of NFI performance. For example, WFA performance is superior to NFI performance for GSS stage 1 patients (Fig. 8)by about 0.20 in ROC AUC (and superior to FFA performance by ∼ 0.06), whereas it is superior to NFI performance across the other stages by ∼0.03). Of note, NFI performed slightly better across the higher GSS stages than FFA, falling intermediate to WFA and FFA at the more advanced stages (Fig. 8) . Thus, it is at the mildest levels of field loss that performance of the methods differs most. 

Finally, performance was considered separately for type of field defect as defined by the GSS to be localized, mixed, or generalized defect. Sixty-one percent of the patients with glaucoma were classified as having only localized defects, 37% had mixed defects, and only one patient had a generalized defect. As seen in Table 4 , whether the defects were mixed or localized had no effect on performance of any of the measures. 

Discussion

In this report, we investigated the ability of quantitative analysis of the shape of the RNFL thickness surface to improve glaucoma detection. The shape analysis technique in general and WFA in particular significantly improved the ability of GDx-VCC polarimetry to detect glaucoma. This improvement over standard VCC methods (NFI) is largest for the earliest cases of glaucoma (stage 1, or MD < −5.0). Specifically, for patients with a stage 1 defect WFA detects 81% of the glaucoma cases at 90% specificity whereas NFI detects only 44% (with FFA intermediate, 61%). 

A prior report from another laboratory ¹¹ applied the Fourier analysis LDF method to data from a VCC prototype device and found sensitivity of 84% at 92% specificity and an ROC area of 0.949 (they also cite an unpublished analysis on a different sample yielding comparable results of 80% sensitivity at 90% specificity and 0.941 ROC area). The study had a similar number of patients and an average MD (−5.91) similar to that in the present study. Those investigators used Fourier analysis amplitude and unsigned phase (and no asymmetry measures). The present results with FFA with the VCC device are quite similar, 0.938 ROC area and 79% sensitivity at 90% specificity using PCA and unsigned phase FFA, and 0.932, 87%/90% with Wilks λ and unsigned phase, which were improved to 0.968 ROC area (90%/90%) and 0.962 (87%/90%), for PCA and Wilks λ, respectively, by adding between-eye asymmetry measures. Use of the WFA method was better than all versions of the FFA method. As reported by Medeiros et al., ¹¹ the performance of Fourier TSNIT shape analysis is generally not as good with GDx data (Essock EA, et al. IOVS 1999;40:ARVO Abstract 3481; Essock EA, et al. IOVS 2001;42:ARVO Abstract 93; Sinai MJ, et al. IOVS 2001;42:ARVO Abstract 717; Sinai MJ, et al. IOVS 2002;43:ARVO E-Abstract 302; Bryant FD, et al. IOVS 2000;41:ARVO Abstract 485; Essock EA, et al. IOVS 2003;44:ARVO E-Abstract 3378) ¹⁷ as with GDx-VCC data, although it is significantly better than standard GDx analysis with the number (i.e., the VCC standard metric, NFI, improves the GDx standard metric, the number, just as TSNIT shape analysis with VCC data improves TSNIT analysis with GDx data). A thorough shape-based analysis with OCT data is yet to be reported, although preliminary analysis is available for an FFA metric. ¹⁷  

Considering LDF analysis of Fourier analysis coefficients alone (i.e., ignoring the superior WFA method), certain conclusions can be drawn about that method from the present and prior results. First, even though the optimal classification varies with the particular sample, a classifier obtained with one sample can be successfully applied to another sample (Essock et al. IOVS 2001;42:ARVO Abstract 93; Bryant et al. IOVS 2000;41:ARVO Abstract 485). Second, certain coefficients (features) tend to appear fairly consistently when trained on different samples from different clinics. For full-curve analysis (i.e., not analyzed by hemifield), amplitudes for the second and fifth harmonics are often the more important features (Essock EA, et al. IOVS 2001;43:ARVO Abstract 93). ^{11 17} In the present study, with a sample of VCC data, the phases for the second, third, and fourth harmonics were most important, after the DC component. The Fourier features that were important for a limited set of OCT data, ¹⁷ DC level, amplitude of the second harmonic, and phase of the third harmonic, fit well with these trends. Clearly what is needed next is to apply FFA, and preferably WFA, to a very large sample of patients with early-stage glaucoma to establish a standardized equation. Another apparent regularity of FFA studies is that performance can be improved significantly by adding interocular asymmetry of RNFL thickness shape (Essock EA, et al. IOVS 1999;40:ARVO Abstract 3481; Bryant FD, et al. IOVS 2000;41:ARVO Abstract 485) ³⁴ (present study; see also Ref. ³⁵ ). All variations of FFA (Wilks λ or PCA, and signed or unsigned phase) had performance improved by adding between-eye asymmetry values. It remains to be seen, however, whether the specific asymmetry values used will remain consistent enough across different samples to be of general use. This study also shows that contrary to prior suggestion (Bryant FD, et al. IOVS 2000;41:ARVO Abstract 485), ³⁴ using within-eye asymmetry of Fourier coefficients does not necessarily improve performance. Finally, whether the signed or unsigned phase is used does not affect performance. This latter observation may suggest that deviations from a TSNIT template matter more than the exact shape of that deviation. WFA, in contrast, is not helped by adding either within-eye or between-eye asymmetry. This suggests that WFA may be more stable across different samples than FFA with between-eye asymmetry measures, in addition to performing the best of all measures tested. 

Due to incorporation of PCA, it is difficult to determine which part of the TSNIT curve, if any, is a critical characteristic in detecting glaucomatous eyes by WFA. Neuroanatomical disruption by glaucoma is of course very idiosyncratic across individuals, and whether any region or shape characteristics are diagnostic of glaucoma in TSNIT analysis is yet to be determined. Perhaps quantifying an eye’s overall deviation from a standard TSNIT “normal template” will be a more fruitful approach. What is apparent here (Fig. 5)is that while the bulk of the difference between glaucomatous and healthy eyes is the well-known diminished amplitude of the superior and inferior humps (and associated decrease of modulation depth), there are other overall changes that are frequent enough to be apparent in the mean curves. Specifically, in the nasal region it appears that a diminution of nasal thickness and a change in the shape of the nasal trough is common. Second, skewing of the superior and inferior peaks, shifting the peaks toward each other (nasally), is also common. That is, of the fibers of the superior and inferior humps, it appears that fibers to the temporal side of these humps are lost more often in glaucoma than fibers to the nasal side of the humps, creating a skewed shape of the superior and inferior humps. Presumably, these are among the shape changes that provide the basis for shape-analysis methods to distinguish glaucomatous eyes from healthy eyes. Two examples of glaucomatous eyes illustrating this nasally shifted skewed shape are shown in Figure 9A and 9B , along with the mean normal data. This feature may prove to be a very useful RNFL shape change and needs to be analyzed in future studies. 

Finally, the jaggedness of the individual TSNIT curves deserves comment (e.g., Figs. 6B 9A 9B ). Pronounced notches of 10 thickness units or more were numerous and were present in about half of the eyes. These notches presumably are not blood vessels as the software employs a blood vessel removal algorithm before producing the TSNIT values. Although some of these could be considered as “split bundles,” these notches occurred at a variety of locations across the TSNIT curves. The presence of notches make it difficult to infer a “typical” shape of the TSNIT curve from group means and also clearly make the job of classifying eyes as healthy or glaucomatous difficult for any shape-based classification method. In the present data, there was a general tendency for jagged TSNIT curves to be misclassified by NFI and FFA more than by WFA. Furthermore, the FFA procedure appeared to be most affected by an asymmetry of height or width between superior and inferior humps. This is consistent with the finding that adding the between-eye asymmetry measure improved FFA performance, whereas WFA performance was not improved by adding this additional measure. These conclusions about curve shape also support the idea that an analysis method like a wavelet analysis that better captures abrupt local change than does a Fourier-based method would be more useful to characterize VCC TSNIT. 

Regarding the limitations of the present study, the primary limitation is that, although the present sample size was large, a very large sample size is needed to obtain a standardized classification function. Relatedly, because the NFI metric was derived previously on a different sample (the company’s normative database), NFI performance would be underestimated in comparison to WFA and FFA performance; the K-fold validation technique used in this study is not effective at preventing the bias. However, simulations performed show that 400 runs are adequate to achieve very stable performance of the validation procedure. Another issue is that, as with most of the literature, focus should be shifted to the very earliest cases of glaucoma, and even preperimetric cases, both of which must be assessed in large numbers. Furthermore, the issue of split bundles per se, and jagged VCC curves in general, should be investigated and quantified. We speculate that adding a WFA parameter to custom “tune” the WFA procedure (or the Fourier based procedure) to individuals based on the presence of multiple superior or inferior humps will help performance considerably. 

Clearly, assessing the classification of patients with early glaucoma is only a first step, and RNFL analysis methods should emphasize the follow-up of patients. Assessment of the risk that field defects may progress in patients with glaucoma, or of the risk that suspected glaucoma may advance and convert to glaucomatous field defects, is of utmost importance. WFA has had initial success assessing such risk using GDx data (Guvant P, et al. IOVS 2004;45:ARVO E-Abstract 5504). ²⁰ This type of work should be extended to VCC data and other RNFL measurement methods with large sample sizes. 

 Supported by a grant from the Kentucky Science and Engineering Foundation (KSEF) and Office of Naval Research Grant 00014-03-1-0224.

 Submitted for publication September 29, 2004; revised February 10 and April 11, 2005; accepted April 14, 2005.

 Disclosure: E.A. Essock, (P); Y. Zheng, (P); P. Gunvant, None

 The publication costs of this article were defrayed in part by page charge payment. This article must therefore be marked “advertisement” in accordance with 18 U.S.C. §1734 solely to indicate this fact.

 

Figure 1.

 

The process of fitting a TSNIT RNFL thickness data set by Fourier analysis is illustrated. Raw data from a healthy eye are shown in the right column. Top to bottom: data set fit by one sine wave (the fundamental), and the sum of the first two, three, and four harmonically related sine waves (F1, F2, F3 and F4, respectively) which are shown in the left column. Optimal sine waves of each spatial frequency (width), phase (relative position), and amplitude are provided and linearly summed, point by point, to approximate the original data.

Figure 2.

 

The scaling function (A) and wavelet function (B) used in the WFA analysis (fourth-order Symlets wavelet).

Figure 3.

 

Flowchart showing the steps involved in classification using the WFA method

Figure 4.

 

WFA performance, defined as the ROC AUC, is plotted as a function of the number of features used.

Table 1.

 

ROC AUC for the Main Analysis Methods

	WFA	FFA (PCA)	FFA (Wilks)	NFI
ROC Area	0.978	0.938	0.932	0.900
Sensitivity at 90	0.925	0.794	0.873	0.716
Sensitivity at 95	0.878	0.729	0.838	0.716

 WFA, wavelet-Fourier analysis; FFA (PCA), fast Fourier analysis with principal component analysis; FFA (Wilks), fast Fourier analysis with Wilks λ; NFI, nerve fiber indicator.

Table 2.

 

ROC AUC with and without Asymmetry Measures

	WFA	FFA (PCA)	FFA (Wilks)	NFI
No asymmetry measures	0.978	0.938	0.932	0.900
Between-eye asymmetry used	0.977	0.968	0.962	0.911
Within-eye asymmetry used	0.954	0.933	0.894	N/A

 Methods are described in Table 1 .

Table 3.

 

Classification Performance of Various Metrics

Dataset	WFA	FFA Wilks λ (Unsigned)	FFA Wilks λ (Signed)	FFA PCA (Unsigned)	FFA PCA (Signed)	NFI
No asymmetry measures	0.892/0.930/0.978 (0.925; 0.878)	0.857/0.908/0.932 (0.873; 0.838)	0.775/0.857/0.909 (0.750; 0.658)	0.835/0.848/0.938 (0.794; 0.729)	0.823/0.842/0.930 (0.779; 0.712)	0.761/0.821/0.900 (0.722; 0.716)
Between-eye asymmetry	0.884/0.946/0.977 (0.927; 0.885)	0.839/0.942/0.962 (0.873; 0.833)	0.895/0.949/0.975 (0.922; 0.894)	0.867/0.969/0.968 (0.903; 0.886)	0.868/0.970/0.969 (0.904; 0.886)	0.717/0.949/0.911† (0.762; 0.709)
Within-eye asymmetry	0.863/0.878/0.954 (0.840; 0.771)	0.761/0.859/0.894 (0.730; 0.630)	0.794/0.861/0.904 (0.761; 0.667)	0.828/0.874/0.933 (0.803; 0.714)	0.828/0.877/0.934 (0.805; 0.715)	N/A

 Signed represents the actual values of phase in the FFA; unsigned represents using the absolute value of phase in FFA. The data report the sensitivity/specificity/ROC area (with ROC AUC in bold type), and the sensitivity is reported in parentheses when specificity is fixed at 90 and at 95. The sensitivity, specificity, and ROC were calculated by averaging performances of 400 runs, except the FFA Wilks, which was calculated by averaging classification performance of 10 runs.

†  This NFI performance was obtained by using a trained LDF with two factors: NFI score and asymmetry (i.e., difference) in NFI score between eyes.

Figure 5.

 

The mean of the reconstructed RNFL patterns from the WFA are shown for healthy (N) and glaucomatous (G) eyes. Mean RNFL patterns reconstructed with backward WFA transform by using only the first (A) or second (B) feature. (C, D) Mean RNFL patterns reconstructed by using the first two (C) or first four (D) features, respectively. The means of raw data from the 67 healthy subjects (solid line) and 67 glaucomatous (dashed line) eyes are also plotted in (D) to facilitate comparison. Note that angles from 0° to 360° indicate the temporal (0°), superior (90°), nasal (180°), inferior (270°), and back to temporal (360°) retina of the eye.

Figure 6.

 

The TSNIT RNFL individual thickness plot (64 sectors) and WFA reconstruction are shown for eyes of two healthy subjects (A, B).

Figure 7.

 

Fit of average data by the two shape-based methods is compared in mean healthy (A) and glaucomatous (B) eyes.

Table 4.

 

Ability of Analysis Methods to Detect Earliest Glaucoma

Grouping Based on Visual Fields	ROC Area
		WFA	FFA	NFI
Glaucoma staging system
Stage 1	0.947	0.888	0.749
Stage 2	0.974	0.921	0.951
Stage 3	0.996	0.960	0.959
Stage 4	0.997	0.952	0.978
Stage 5	0.983	0.983	0.949
Grouped by MD
Group 0 to −5	0.961	0.903	0.816
Group −5 to −10	0.990	0.958	0.964
Group −10 to −15	0.997	0.956	0.976
Group −15 to −20	0.979	0.982	0.934
Group −20 to −25	1.000	0.946	1.000
Type of defect
Localized	0.977	0.941	0.892
Mixed	0.980	0.946	0.909
Generalized	0.983	0.536	1.000

Figure 8.

 

The performances of the three metrics (WFA, NFA, and NFI) are illustrated as a function of severity of field loss (GSS stage).

Figure 9.

 

Two examples (A, B) of RNFL thickness of glaucomatous eyes are shown. In comparison to the mean normal curve, the superior and inferior peaks are shifted toward each other.

The authors thank Bill Rising for statistical consultations.