April 2007
Volume 48, Issue 4
Free
Glaucoma  |   April 2007
Automated Segmentation of the Optic Disc from Stereo Color Photographs Using Physiologically Plausible Features
Author Affiliations
  • Michael D. Abràmoff
    From the Department of Ophthalmology and Visual Sciences, University of Iowa Hospitals and Clinics, Iowa City, Iowa; the
    Department of Veterans Affairs, Iowa City VA Medical Center, Iowa City, Iowa; and the
    Department of Electrical and Computer Engineering, University of Iowa, Iowa City, Iowa.
  • Wallace L. M. Alward
    From the Department of Ophthalmology and Visual Sciences, University of Iowa Hospitals and Clinics, Iowa City, Iowa; the
    Department of Veterans Affairs, Iowa City VA Medical Center, Iowa City, Iowa; and the
  • Emily C. Greenlee
    From the Department of Ophthalmology and Visual Sciences, University of Iowa Hospitals and Clinics, Iowa City, Iowa; the
    Department of Veterans Affairs, Iowa City VA Medical Center, Iowa City, Iowa; and the
  • Lesya Shuba
    From the Department of Ophthalmology and Visual Sciences, University of Iowa Hospitals and Clinics, Iowa City, Iowa; the
  • Chan Y. Kim
    From the Department of Ophthalmology and Visual Sciences, University of Iowa Hospitals and Clinics, Iowa City, Iowa; the
  • John H. Fingert
    From the Department of Ophthalmology and Visual Sciences, University of Iowa Hospitals and Clinics, Iowa City, Iowa; the
  • Young H. Kwon
    From the Department of Ophthalmology and Visual Sciences, University of Iowa Hospitals and Clinics, Iowa City, Iowa; the
Investigative Ophthalmology & Visual Science April 2007, Vol.48, 1665-1673. doi:10.1167/iovs.06-1081
  • Views
  • PDF
  • Share
  • Tools
    • Alerts
      ×
      This feature is available to Subscribers Only
      Sign In or Create an Account ×
    • Get Citation

      Michael D. Abràmoff, Wallace L. M. Alward, Emily C. Greenlee, Lesya Shuba, Chan Y. Kim, John H. Fingert, Young H. Kwon; Automated Segmentation of the Optic Disc from Stereo Color Photographs Using Physiologically Plausible Features. Invest. Ophthalmol. Vis. Sci. 2007;48(4):1665-1673. doi: 10.1167/iovs.06-1081.

      Download citation file:


      © 2015 Association for Research in Vision and Ophthalmology.

      ×
  • Supplements
Abstract

purpose. To evaluate a novel automated segmentation algorithm for cup-to-disc segmentation from stereo color photographs of patients with glaucoma for the measurement of glaucoma progression.

methods. Stereo color photographs of the optic disc were obtained by using a fixed stereo-base fundus camera in 58 eyes of 58 patients with suspected or open-angle glaucoma. Manual planimetry was performed by three glaucoma faculty members to delineate a reference standard rim and cup segmentation of all stereo pairs and by three glaucoma fellows as well. Pixel feature classification was evaluated on the stereo pairs and corresponding reference standard, by using feature computation based on simulation of photoreceptor color opponency and visual cortex simple and complex cells. An optimal subset of 12 features was used to segment all pixels in all stereo pairs, and the percentage of pixels assigned the correct class and linear cup-to-disc ratio (LCDR) estimates of the glaucoma fellows and the algorithm were compared to the reference standard.

results. The algorithm was able to assign cup, rim, and background correctly to 88% of all pixels. Correlations of the LCDR estimates of glaucoma fellows with those of the reference standard were 0.73 (95% CI, 0.58–0.83), 0.81 (95% CI, 0.70–0.89), and 0.86 (95% CI, 0.78–0.91), respectively, whereas the correlation of the algorithm with the reference standard was 0.93 (95% CI, 0.89–0.96; n = 58).

conclusions. The pixel feature classification algorithm allows objective segmentation of the optic disc from conventional color stereo photographs automatically without human input. The performance of the disc segmentation and LCDR calculation of the algorithm was comparable to that of glaucoma fellows in training and is promising for objective evaluation of optic disc cupping.

Glaucoma is the second leading cause of blindness in the United States, characterized by gradual damage of the optic nerve and resultant visual field loss. 1 Early diagnosis and optimal treatment have been shown to minimize the risk of visual loss due to glaucoma. 2 The hallmark of glaucomatous progression is cupping of the optic disc. One way of determining the amount of cupping is planimetry performed by experienced glaucoma specialists examining stereo color photographs of the optic disc. 3 4  
Pixel feature classification is a machine learning technique that can be used to assign a class to the pixels in an image, in this case the cup, rim, or background. 5 6 7 8 Pixel feature classification uses multiple pixel features, which are numeric properties of a pixel and its surroundings. Originally, pixel intensity was used, and later, its contrast with the surrounding region, its proximity to an edge, and so forth, were added. As a supervised algorithm, it has two stages: a training stage, in which the algorithm “statistically learns” to classify pixels correctly from known classifications, and a testing or classification stage, in which the algorithm classifies a new, previously unseen image. The number of features for any pixel is essentially infinite, and it is easy to prove that one or more subsets of this infinite set are optimal for classifying the image according to some reference standard. Hundreds of features for a pixel can be calculated in the training stage to cast as wide a net as possible. Still, because all potential features cannot possibly be taken into account, such a selection is arbitrary to some extent and carries a bias. One of the advantages of the pixel feature classification method is that it avoids the potentially much larger bias of an imaging expert’s developing the algorithm by selecting one or two feature detectors. This is common in unsupervised image classification, such as in an edge-based, deformable model or in template-based algorithms that have been used in segmentation of retinal images and optic disc in nonstereo images. 9 Pixel feature classification has so far outperformed all other approaches for the segmentation of vessels, hemorrhages and microaneurysms, exudates, infarcts, and drusen in retinal images 7 8 (Niemeijer et al. IOVS 2005;45:ARVO E-Abstract 3468). 
Another advantage of pixel feature classification is its potential for mimicking human experts’ visual performance by using parallel visual “channels.” 10 We have designed pixel features for color that mimic the output of the dark-bright, red-green, and blue-yellow center surround opponency cells, with different receptive field sizes and orientation, and features for stereo disparity that mimic the output of binocular cells in area V1 of the visual cortex. 11 12  
The purpose of this study was to compare the performance of pixel feature classification that segments the optic disc into cup, rim, and background from color stereo photographs and the performance of glaucoma fellows with the results of three glaucoma faculty experts, for measurement of glaucoma progression. 
Specifically, we compared the performance of the pixel classification algorithm to segment the optic disc photographs into cup, rim, and background with the planimetric performance of the fellows against the reference standard planimetric segmentation by the glaucoma experts. 
Methods
Stereo Color Images of the Disc
The study was approved by the University of Iowa’s Institutional Review Board and adhered to the tenets of the Declaration of Helsinki. For the test set, stereo photographs of the optic disc were obtained from 58 consecutive patients of the Glaucoma Service at the university. Patients were included if they met the definition of suspected glaucoma or open-angle glaucoma. Those with a history of angle closure or combined-mechanism glaucoma, or any nonglaucomatous optic neuropathy were excluded. The diagnosis of “suspected glaucoma” was based on a suspect optic nerve appearance (enlarged cupping on clinical examination) with normal visual field and intraocular pressure (IOP) ≤21 mm Hg. It also included those with normal optic disc appearance on biomicroscopy and normal visual field, but with elevated IOP (>21 mm Hg). The “open-angle glaucoma” diagnosis included patients with both primary and secondary open-angle glaucoma and was based on the presence of an open iridocorneal angle, glaucomatous optic disc, and/or nerve fiber layer defects on biomicroscopy, and visual field changes (regardless of IOP level). Glaucomatous optic discs were identified as those with either diffuse or focal thinning of the neuroretinal rim. Visual field abnormalities were considered to be glaucomatous if they were consistent with the optic nerve examination and had either a (1) typical nerve fiber layer distribution or (2) glaucoma hemifield test result outside normal limits. The diagnoses were made by the treating faculty glaucoma specialist according to these definitions. 
Color slide stereo photographs centered on the optic disc of both eyes were acquired with a fixed stereo-base retinal camera (3Dx; Nidek, Newark, NJ). After development, slides were scanned at a 4096 × 4096-pixel resolution, 24-bit depth with a slide scanner (Kodak Inc., Rochester, NY). Only one eye of each patient was selected for the study based on the quality of the stereo image. To limit the number of pixels (>16 million for each left–right stereo pair) that would otherwise have to be processed by the algorithm, the left and right stereo images of the selected stereo pairs were cropped to 512 × 512 pixels each with the optic disc in the center, by automatically locating it in the 4096 × 4096 images and cropping to 256 pixels outward from the center of the disc based on its brightness. 13 The cropped images included the optic disc in its entirety in all images (see Fig. 1 ). 
Computer-aided planimetry of the stereo pairs with a hand-held stereoscope was performed by three faculty glaucoma specialists at the University of Iowa on the 58 stereo pairs by using custom-made annotation software. 14 In addition, three glaucoma fellows in the university’s Glaucoma Fellowship performed planimetry on all stereo pairs in the same manner. The human graders were masked from the patient identifier, and the clinical information was associated with the photographs and the results of all other graders at the time of grading. Graders were instructed to outline the disc margin while viewing the photographs with a stereo viewer. Then, they were instructed to draw the cup to the best of their ability. When the cup margin was difficult to identify (usually at the temporal slope of the disc), they were instructed to identify and draw the cup margin on the slope that is half way between the bottom of the cup and top of the rim. 
A reference standard grading with a unique class (cup, rim, and background) assigned to each pixel was then created by combining the grading by each faculty grader in a “majority-win” manner: each pixel was assigned the class that received the majority of votes (e.g., if two graders voted a pixel to be rim and one voted it to be cup, the pixel assigned was rim). In the case of a draw, the pixel was assigned the background class (Figure 2)
Pixel Feature Classification
Rationale for Inclusion of Specific Features.
Originally, pixel intensity and simple edge operators were used as pixel features. 5 Once it was realized that a pixel’s feature could also incorporate the adjacent pixels (i.e., its surround) and the Gaussian filter bank concept had been established mathematically, 15 it was possible to use filter bank outputs to create a larger set of features, from which an optimal combination of features could be selected—a method used initially in pulmonary image analysis. 16 Gaussian filter bank features are sensitive to edges and textures in the gray-scale intensity image obtained from a color image at different scales and orientations. For more optimal use of the color information in fundus color images, we added the outputs of a Gaussian steerable (i.e., having multiple orientations) filter bank in hue, saturation, brightness space as features and also the variance in the red, green, and blue channels in a small region around the pixel. 8  
Such features are obviously different from the color opponency model of primate color vision, and for this study we added the output of a Gaussian steerable filter bank in color opponency space for dark-bright, red-green, and blue-yellow center-surround color opponency. See Appendix A for the mathematical details. Preliminary experiments showed that the color opponency filter bank improved the performance (Figure 3)
Similarly, because expert human graders require stereo photographs for optimal segmentation of the optic nerve head, we decided to add horizontal stereo disparity maps as features. In the literature, many known algorithms use edge and texture information and some form of energy minimization to solve this depth from stereo problem in color or video images. 17 These algorithms are fundamentally different from what is known about stereo disparity computation in primates. 18 In preliminary studies, we found that the stereo disparity algorithms that have been found to be optimal for natural scene stereo disparity—namely, the graph search algorithms first described by Kolmogorov and Zabih 17 and Birchfield and Tomasi 19 —did not perform as well as expected on the low-contrast/high-noise images in this study. Therefore, we decided to add features calculated by simulating the output of binocular complex cells in V1 visual cortex as first described by Qian and Zhu. 12 The output of binocular simple cells in the mammalian V1 has been shown to correspond closely with the convolution of its receptive field with Gabor-basis kernels, with the same spatial frequency and scale but different phases for the left and right fields. Stereo disparity–sensitive complex cells can then be modeled by summing the squared outputs of two simple cells with different phases, but the same size receptive field. See Appendix B for the mathematical background. 
An informal taxonomy of pixel features obtainable from stereo color images may be:
  •  
    “Less physiologically plausible”: The outputs of edge detectors such as those described by Canny and Sobel, or graph search stereo disparity algorithms, as there is currently no evidence that these are expressed by any neural circuit. 20
  •  
    “Somewhat physiologically plausible”: The outputs of Gaussian filter banks in intensity, or RGB, space, because there is evidence that a Gaussian profile is expressed by certain neural circuits. 21 No red/green/blue– or hue/saturation/intensity–sensitive receptive fields have been described to our knowledge.
  •  
    “More physiologically plausible”: Gabor wavelets in color opponency space, and the output of simulations of binocular complex cells, as there is evidence that these are expressed in certain neural circuits. 22 23
We have assumed that a combination of retinal and cortical image processing as found in biological systems would be effective in classifying the stereo color images of the optic disc and thus have implemented these physiologically plausible features and added them to the set of potential features from which an optimal set of features is then selected. Whether these features indeed increase the performance of the system can be studied (implicitly) by running feature selection on sets with and without the “more physiologically plausible” features and comparing the accuracy. 
A Priori Expectation Feature
In pixel feature classification, non pixel-based information is not naturally included. However, it is easy to see that in stereo disc color imaging of this type, the likelihood that a pixel in the center of the image is cup is higher than the likelihood that it is rim or background. To capture this a priori expectation of the class of a pixel based on its location in the image, an a priori rim–nonrim feature was calculated by performing principal component analysis on all reference standard images and taking the first principal component. 24 This feature expresses the expectation of any pixel’s being rim or either background or cup based on its position. Preliminary studies without the a priori feature did not attain an accuracy of more than 73% (see the Results section). 
A total of 108 + 3 × 48 + 1 features were thus calculated, for all 15,204,352 pixels in all 58 color stereo pairs. Because these features can be calculated independently in parallel, as is the case in the primate visual system, this is a naturally parallel computation. See Figure 4for an example of all values of the 12 optimal pixel features that were selected for a specific stereo pair. 
Classifier Implementation
In pixel feature classification, a classifier is trained on a subset of all data, the training set, and then repeatedly tested on another subset of the data, the test set. By repeatedly testing with different subsets of all the features and comparing the performance of the classifier with these different feature subsets, the optimal subset of features can be selected. In preliminary studies, the k-nearest neighbor (kNN) classifier outperformed support vector machines and linear discriminant analyzer classifiers, as has been found in other studies. 8 16 25 26 A kNN classifier has the advantage that no elaborate training is required and that only one parameter needs to be adjusted, the k. We have interfaced the publicly available kNN library (www.cs.umd.edu/∼mount/ANN/ developed by David Mount and Sunil Arya and provided in the public domain by the Department of Computer Science, University of Maryland, College, Park, MD) to Java. 27  
Data Analysis
The classifier was optimized on the accuracy as explained in Appendix C to find the set of optimally discriminant features. A clinically more relevant measure was also used to compare the performance of the pixel feature classification and the three glaucoma fellows to the reference standard. For this purpose, all pixels in all stereo pairs were classified in a leave-one-out fashion. 25 The resultant gradings contained the algorithm’s estimate for each pixel as to whether that pixel was cup, rim, or background, similar to the gradings by the glaucoma faculty and the glaucoma fellows. A linear cup-to-disc ratio was then computed for each graded image as a square root of area cup-to-disc ratio:  
\[cd_{\mathrm{r,linear}}{=}\sqrt{\frac{n_{\mathrm{c}}}{n_{\mathrm{r}}{+}n_{\mathrm{c}}}},\]
where n c is the number of pixels in cup, and n r is the number of pixels in rim. 28  
Optimal Feature Selection
The algorithm had been trained on a larger set of 101 stereo disc pairs of patients with glaucoma (Abràmoff et al., IOVS 2006;46:ARVO E-Abstract 3635). The training set was used to determine which combination of features gives the best accuracy (“feature selection,” see Appendix C). A hard classification using “majority-win” of the kNNs was used as the classification for each pixel. Each feature column in the training set was scaled to have zero mean and unit variance over all pixels. The most discriminant features were then selected in a sequential forward-floating search, 25 by adding all features to an initially empty feature set one by one and retaining only the new feature that showed the highest gain in classification performance. Preliminary studies showed that k = 501 gave slightly better performance over a range of feature vector lengths, but sensitivity to changing from 1 to 1001 was low. The number of samples was 1200 per image, per class. The pixel feature classification was implemented in Java, and it took 15 minutes to calculate 12 features and completely grade a stereo pair on a standard 2-GHz 2-Gbyte computer (running Windows XP; Microsoft, Redmond, WA). 
Table 1shows the 12 features that were found to be most dominant. Figure 5shows an example stereo pair and those pixel features that were most discriminant features. After feature selection, the algorithm was ready for testing using this optimal set of features. 
Results
The stereo pairs of 58 eyes of 58 patients were included to test the algorithm. The age of these patients was 61.2 ± 17.8 years (mean ± SD), and there were 33 women. There were 22 patients with suspected glaucoma and 36 with open-angle glaucoma. Of the total of 58 eyes, 27 underwent Humphrey visual field testing (24-2 SITA Standard, HFA II; Carl Zeiss Meditec, Dublin, CA) within 1 year of the disc photographs. The remainder (n = 31) had either Goldmann visual field testing (n = 28, Goldmann Perimeter 940; Haag Streit, Liebefeld, Switzerland) or no field tests (n = 3). We performed Goldmann visual field tests in patients whose Humphrey visual field result was not reliable and/or who had moderate to severe visual field defect (typically −15 dB mean deviation [MD] or worse). In these patients, we found the Goldmann visual field results to be more reliable. A few patients (n = 3) were unable to perform any visual field tests at all (e.g., those with Alzheimer’s disease). The Humphrey visual field mean defect (in decibels ± SD) in eyes with suspected glaucoma (n = 11) or open-angle glaucoma (n = 16) were −1.1 ± 2.3 and −4.1 ± 5.2, respectively. The Humphrey visual field pattern SD (dB ± SD) in eyes with suspected glaucoma and open-angle glaucoma was 2.8 ± 3.6 and 4.5 ± 4.5, respectively. The Goldmann visual field results of the suspected-glaucoma group (n = 8) were within normal limits, but the results in the open-angle group (n = 20) showed a variety of defects consistent with glaucomatous nerve fiber layer loss including paracentral scotoma, nasal step, and arcuate defect(s). 
Average disc size (cup and rim area, as determined from the reference standard) was 14,085 ± 3,467 (SD) pixels. All pixels in all stereo pairs were graded by the algorithm. Accuracy was 0.876; the algorithm was able to assign 88% of all pixels in the 58 stereo pairs to the correct class. Figure 5shows how the gradings by the glaucoma faculty, the glaucoma fellows, and the algorithm compared with each other visually. Correlations of the linear cup to disc ratio for each of the glaucoma fellows with the reference standard were 0.73 (95% confidence interval [CI], 0.58–0.83), 0.81 (95% CI, 0.70–0.89), and 0.86 (95% CI, 0.78–0.91), respectively, whereas correlation of the algorithm with the reference standard was 0.93 (95% CI, 0.89–0.96; n = 58). 
Subgroup analysis showed no significant differences in the correlation between the suspected and open-angle glaucoma groups for each of the fellows or the algorithm, and there was no dependency on disc size (P > 0.05). Figure 6shows scatterplots of the estimates for the linear cup-to-disc ratio (CDR) for each stereo pair for the glaucoma fellows and the algorithm versus the reference standard. 
Discussion
Our results show that pixel feature classification with biologically inspired features can classify a set of optic disc stereo images with an accuracy that is comparable to that of glaucoma fellows. Other studies have also shown agreement of pixel classification algorithms for the segmentation of retinal vessels and the detection of diabetic retinopathy lesions with classifications made by human experts. 8 26 In fact, Niemeijer et al., 7 found that a pixel classification–based approach for retinal vessel segmentation outperformed other supervised and unsupervised methods. The same group of investigators used a pixel classification approach to detect red diabetic retinopathy lesions and reached a sensitivity of 100%, a specificity of 87%, and an area under the ROC curve of 98.5% in a study of 100 patients. 8  
We are aware of only one previous study of classification of the optic disc in which stereo photographs were used. 29 In that study, a matching cost approach was used, based on pyramidal cross correlation. Recent depth from stereo algorithms, mostly based on graph searches, have been shown to outperform matching cost computation algorithms on real-world depth scenes. 17 19 30 As shown by our results, these graph search algorithms in their turn can be outperformed by an algorithm that mimics binocular cells in the visual cortex in optic disc stereo images. 12  
On a computational level, the role of most features is to determine how much a patch around the pixel resembles the filter kernel of that feature. The results indicate that the features selected for optimal accuracy of the system were predominantly of the “somewhat” and “more” physiologically plausible categories, and not for example, the familiar Sobel or Canny intensity edge detectors. On an image-understanding level, it is difficult to explain the role of the selected features. A second-order Gaussian derivative of a specific orientation and blue-yellow color opponency can be said to correspond to “striped oriented patches” of the pixel surround, but is more than a simple edge detector, and also different from a texture detector. This situation is analogous to the difficulty in describing the role of different receptive field properties for scene analysis by the primate visual cortex. 
Objective classification of optic disc cup and rim from color stereo photographs is promising, because it allows optic disc stereo photographs to be evaluated objectively and rapidly. However, the ultimate goal of an algorithm such as this is the objective measurement of glaucoma progression, or enlargement of the optic cup, over time. The results of the present study are intended foremost to serve as a proof of concept, and many aspects can be improved, including adding features that mimic physiology even more and adding vessel segmentation. 
Nevertheless, several serious issues remain to be addressed before the present algorithm can be thus used. First, the algorithm is deterministic and therefore returns the same classification if offered the same stereo pair repeatedly. However, whether it also gives very similar classifications on two stereo pairs of the same patient of the same eye taken minutes apart, with subsequent patient movement, vignetting, and refocusing, has not been determined. 
Second, there is room for substantial improvement, as is evident from the noisy visual appearance of the classification, compared with the smooth boundaries of the reference standard, and the percentage of correctly classified pixels, which is only 88%. For example, compared with the optic disc images obtained with a fundus camera (Carl Zeiss Meditec), the 3Dx (Nidek) gives decidedly darker and less crisp images, which may have affected the quality of the performance negatively. We have adopted several different postprocessing techniques that are showing promising results. 31  
Third, the algorithm performs well on stereo pairs obtained with the fixed-base 3Dx, which may not be widely available, and it has not been established whether the algorithm will perform similarly on non–fixed-base stereo pair cameras. 
Fourth, the algorithm’s performance is highly dependent on a single feature, the prior PCA. One of its roles is likely to be to indicate that the disc is in the center of the image, and the cup is in the center of the disc and surrounded by the rim. Because this is an assumption about the expected size and location of the disc, the algorithm’s accuracy may be lower on discs that are smaller or larger than in our sample and also if the disc is not properly centered in the cropped image. Therefore, we are currently engaged in removing this feature from the selected features while maintaining performance. 
Finally, the algorithm is slow. It is important to understand that the current implementation of the algorithm was not written for speed, but for ease of use and modifiability, and that it can easily be optimized for speed. In our previous optimized algorithms, classification times of 15 seconds for a whole image were reached. 8 In addition, calculation of the features and feature classification is naturally parallel, so that with appropriate hardware, the features for all pixels can be calculated in parallel. 
Contrary to expectations, stereo disparity was not the most discriminant feature for the performance of the algorithm. In Figure 1 , it can be appreciated that the stereo photographs are not very crisp. Possibly, using a pair of photographs improves the quality of the fused image by suppressing these irregularities to some extent. In addition, human graders may be unaware of the relative importance of color and texture of the tissues compared with depth-from-stereo as the feature they use to segment the optic disc. The algorithm may therefore also perform reasonably on non-stereo color photographs. 
In summary, the results of this study have shown that pixel feature classification with biologically inspired features of stereo color images is a good starting point for the classification of the optic disc and other three-dimensional structures on the retina, and also, that features based on the physiological process outperform standard pixel features on segmentation of the optic disc. 
Appendix 1
Color Opponency Gaussian Filter Bank
By careful selection of steerable Gaussian kernels in a Gaussian filter bank, second-order invariant features can be calculated, so that all other features up to the second order are combinations of the features in this invariant set. 15 Color opponency Gaussian filter banks 32 L db,L rg,L by, for dark-bright (db), red-green (rg), and blue-yellow (by) opponency images, are computed from each color image as follows:  
\[L_{\mathrm{db}}{=}\ \frac{L_{\mathrm{r}}{+}L_{\mathrm{g}}{+}L_{\mathrm{b}}}{3},\]
 
\[L_{\mathrm{rg}}{=}L_{\mathrm{r}}{-}L_{\mathrm{g}},\]
and  
\[L_{\mathrm{by}}{=}(L_{\mathrm{r}}{+}L_{\mathrm{g}}){-}2L_{\mathrm{b}},\]
with L r,L g,L b the red (r), green (g), and blue (b) channels of the left (L) image of the stereo pair, respectively. 
L db,L rg,L by are then each convolved with Gaussian filter bank kernels to obtain the following color-opponency Gaussian derivative features:  
\[L_{0}({\sigma}),L_{1}^{0{^\circ}}({\sigma}),L_{1}^{90{^\circ}}({\sigma}),L_{2}^{0{^\circ}}({\sigma}),L_{2}^{60{^\circ}}({\sigma}),L_{2}^{120{^\circ}}({\sigma}),\]
with L n,op α(σ) the feature obtained by computing L n α(σ)=L 1G n α(σ) for order n, orientation α (0° and 90° orientation to the local gradient for first-order derivatives and 0°, 60°, and 120° to the local gradient for second-order derivatives), scale σ∈(2,4,8,16,32,64), and opponency (op) image L op∈ (L db,L rg,L by). ⊗ denotes convolution. L n α is computed as follows:  
\[\begin{array}{l}L_{0}{=}L_{1}\\L_{1}^{{\alpha}}{=}\mathrm{cos}{\alpha}L_{x}{+}\mathrm{sin}{\alpha}L_{x}\\L_{2}^{{\alpha}}{=}\mathrm{cos}^{2}{\alpha}L_{xx}{+}\mathrm{cos}{\alpha}\mathrm{sin}{\alpha}L_{xy}{+}\mathrm{sin}^{2}{\alpha}L_{yy},\end{array}\]
where the xx, yy, and xy subscripts to L denote differentiation of L op to x and/or y 15 and results in 108 color opponency Gaussian filter bank features. 
Appendix 2
Stereo Features
Gabor features G σ,op α were calculated by convolution of the three opponency images op∈(L db,L rg,L by) with Gabor wavelet kernels at scales σ = 1, 2, 4, 8, 16, and 32 pixels, at orientation α = 0°, 45°, 90°, and 135°:  
\[\mathit{G}_{{\sigma},op}^{{\alpha}}(x,y){=}e^{{-}\ \frac{x^{2}}{2{\sigma}^{2}}}cos\ ({\omega}x{+}{\varphi}),\]
with phase φ 0 and with spatial frequency ω=
\(\frac{{\pi}}{{\sigma}}\)
constant, for a total of 82 features. See Figure 3for an example of these kernels. Stereo disparity features were computed by first convolving the left and right color opponency images with left and right Gabor kernels k l and k r:  
\[k_{l}(x){=}e^{{-}\ \frac{x^{2}}{2{\sigma}^{2}}}\mathrm{cos}({\omega}x{+}{\varphi}_{\mathrm{l}})\]
 
\[k_{\mathrm{r}}(x){=}e^{{-}\ \frac{x^{2}}{2{\sigma}^{2}}}\mathrm{cos}({\omega}x{+}{\varphi}_{\mathrm{r}}),\]
where φlr are the phases of the left and right Gabor wavelets, and ω its preferred spatial frequency, which was kept constant at ω=  
\[\frac{{\pi}}{{\sigma}}\]
. The linear combination of convolutions of the left and right images with these respective kernels mimics the output of a simple cell:  
\[S_{{\omega}}({\sigma},{\varphi}_{\mathrm{l}},{\varphi}_{\mathrm{r}}){=}L_{\mathrm{l,1}}{\otimes}k_{\mathrm{l}}{+}L_{\mathrm{l,r}}{\otimes}k_{\mathrm{r}},\]
with L 1,l,L l,r the left and right intensity images, respectively. The sum of the squares of two  
\[S_{1,{\omega}}({\sigma},{\varphi}_{\mathrm{l}},{\varphi}_{\mathrm{r}}),S_{2,{\omega}}({\sigma},{\varphi}_{1}{+}\ \frac{{\pi}}{2},{\varphi}_{\mathrm{r}}{+}\ \frac{{\pi}}{2})\]
formed quadrature pairs:  
\[C_{{\omega}}({\sigma},{\Delta}{\varphi}){=}(S_{\mathrm{1,{\omega}}}({\sigma},{\varphi}_{\mathrm{l}},{\varphi}_{\mathrm{r}}))^{2}{+}\left[S_{2,{\omega}}\left({\sigma},{\varphi}_{\mathrm{l}}{+}\ \frac{{\pi}}{2},{\varphi}_{\mathrm{r}}{+}\ \frac{{\pi}}{2}\right)\right]\ ^{2},\]
with Δφ the phase difference  
\[{\Delta}{\varphi}{=}{\varphi}_{\mathrm{r}}{-}{\varphi}_{\mathrm{l}}\]
between the right and left eye receptive fields. With  
\[{\Delta}{\varphi}{\in}\left\{{\pi}\ \frac{6}{8}{\pi}\ \frac{4}{8}{\pi},\ \frac{2}{8}{\pi},0,{-}\ \frac{2}{8}{\pi},{-}\ \frac{4}{8}{\pi},{-}\ \frac{6}{8}{\pi}\right\}\]
, C ω(σ,Δφ) was then tuned to stereo disparities of {−4,−3,−2,−1,0,1,2,3} pixels, respectively, calculated at scales σ =1, 2, 4, 8, 16, and 32 pixels, and for each of the (L op,R op)∈{(L db,R db),(L rg,R rg),(L by,R by)} color opponency pairs, resulting in 48 features for each of the color opponency spaces. 12  
Appendix 3
Balanced Accuracy
Performance was calculated in the form of balanced accuracy (a combination of sensitivity and specificity in a single number) by calculating the number of correctly classified pixels in equally sized samples from the rim, the cup, and the background  
\[{\pi}{=}\ \frac{N_{\mathrm{r,T}}{+}N_{\mathrm{c,T}}{+}N_{\mathrm{b,T}}}{N_{\mathrm{r}}{+}N_{\mathrm{c}}{+}N_{\mathrm{b}}},\]
with N r,T the number of pixels correctly classified in a sample of rim pixels, N c,T the number of pixels classified correctly in a sample of cup pixels, N b,T the number pixels classified correctly as background in a sample of cup pixels, and N r=N c=N b the total number of pixels in each sample. The reason to make this calculation is that otherwise, because of the predominance of background pixels in the image, an algorithm that is “good” in classifying background but “bad” in classifying rim would have an unfair advantage. 
 
Figure 1.
 
Stereo color image pair of the optic nerve head of the left eye after alignment and cropping to 512 × 512 pixels. The left and right stereo images can be fused by holding the page at approximately 30 cm (12 in.) distance from the eye and gazing into infinity.
Figure 1.
 
Stereo color image pair of the optic nerve head of the left eye after alignment and cropping to 512 × 512 pixels. The left and right stereo images can be fused by holding the page at approximately 30 cm (12 in.) distance from the eye and gazing into infinity.
Figure 2.
 
(AC) Gradings of an optic disc stereo pair by three faculty glaucoma specialists. Rim is grayish, cup is whitish in the left image of the stereo pair from Figure 1 . (D) Reference standard created from (AC) with white, cup; gray, rim; and black, background
Figure 2.
 
(AC) Gradings of an optic disc stereo pair by three faculty glaucoma specialists. Rim is grayish, cup is whitish in the left image of the stereo pair from Figure 1 . (D) Reference standard created from (AC) with white, cup; gray, rim; and black, background
Figure 3.
 
Example of the variety of color opponency, steerable, Gaussian filter bank kernels. Kernels are oriented as if the local gradient were horizontal. First row, from left to right: dark-bright opponency kernels for zeroth order, first order 0° to the local gradient, first-order 90° to the local gradient, second-order 0° to the local gradient, second order 60° to the local gradient, and second order 150° to the local gradient, at scale σ 32 pixels. Second row: same, for scale σ 64 pixels; third row: 128 pixels. Next three rows: same, for blue-yellow opponency kernels; last three rows: same, for red-green kernels. Smaller scales are not shown because they are difficult to see. These kernel images also represent the response of each of the features to an impulse function.
Figure 3.
 
Example of the variety of color opponency, steerable, Gaussian filter bank kernels. Kernels are oriented as if the local gradient were horizontal. First row, from left to right: dark-bright opponency kernels for zeroth order, first order 0° to the local gradient, first-order 90° to the local gradient, second-order 0° to the local gradient, second order 60° to the local gradient, and second order 150° to the local gradient, at scale σ 32 pixels. Second row: same, for scale σ 64 pixels; third row: 128 pixels. Next three rows: same, for blue-yellow opponency kernels; last three rows: same, for red-green kernels. Smaller scales are not shown because they are difficult to see. These kernel images also represent the response of each of the features to an impulse function.
Figure 4.
 
Stereo pair and its most discriminant pixel features. (A, B) The left and right stereo image of the stereo pair. Rows below them show images of the most dominant 12 features in rank order (see Table 1for a description of these features).
Figure 4.
 
Stereo pair and its most discriminant pixel features. (A, B) The left and right stereo image of the stereo pair. Rows below them show images of the most dominant 12 features in rank order (see Table 1for a description of these features).
Table 1.
 
Twelve Features Selected as the Most Dominant
Table 1.
 
Twelve Features Selected as the Most Dominant
Rank Description Scale (Pixels) Cumulative Accuracy on Training Set
1 A priori probability of pixel being cup 1 0.810
2 Red-green opponency (opp.) second-order Gaussian derivative transversal to local gradient 32 0.825
3 A priori probability of a pixel being rim 1 0.829
4 Red-green opp. Gaussian blurred (zeroth order) 8 0.832
5 Complex cell output tuned to a disparity of −12 pixels 16 0.845
6 Blue-yellow opp. second-order Gaussian derivative 32 0.854
7 Red-green opp. first-order Gaussian derivative perpendicular to local gradient 16 0.869
8 Bright dark opp. Gaussian blurred (zeroth order) 16 0.872
9 Red-green opp. first-order Gaussian derivative radial to local gradient 4 0.874
10 Blue-yellow Gaussian blurred (first-order) 4 0.874
11 Red-green opp. first-order Gaussian derivative parallel to local gradient 4 0.875
12 Red-green opp. Gaussian blurred (zeroth order) 32 0.876
Figure 5.
 
Classification of stereo pairs by glaucoma faculty, three glaucoma fellows, and pixel feature classification. Rows 1–4: small, medium, and large disc excavation and excavation with inferior notching, respectively. From left to right: left and right images of the stereo pair; reference standard by three glaucoma expert faculty; grading by fellow A, fellow B, and fellow C; and grading by pixel feature classification.
Figure 5.
 
Classification of stereo pairs by glaucoma faculty, three glaucoma fellows, and pixel feature classification. Rows 1–4: small, medium, and large disc excavation and excavation with inferior notching, respectively. From left to right: left and right images of the stereo pair; reference standard by three glaucoma expert faculty; grading by fellow A, fellow B, and fellow C; and grading by pixel feature classification.
Figure 6.
 
Scatterplots of linear cup-to-disc ratio estimates from all 58 stereo pairs by three glaucoma fellows and the pixel feature classification algorithm (vertical axes), against the reference standard grading by three glaucoma faculty (horizontal axes). Clockwise, fellows A and B, the algorithm, and fellow C respectively. r, is the correlation coefficient for each of these estimates with the reference standard.
Figure 6.
 
Scatterplots of linear cup-to-disc ratio estimates from all 58 stereo pairs by three glaucoma fellows and the pixel feature classification algorithm (vertical axes), against the reference standard grading by three glaucoma faculty (horizontal axes). Clockwise, fellows A and B, the algorithm, and fellow C respectively. r, is the correlation coefficient for each of these estimates with the reference standard.
TielschJM, SommerA, KatzJ, RoyallRM, QuigleyHA, JavittJ. Racial variations in the prevalence of primary open-angle glaucoma. The Baltimore Eye Survey. JAMA. 1991;266:369–374. [CrossRef] [PubMed]
HeijlA, LeskeMC, BengtssonB, HymanL, BengtssonB, HusseinM. Reduction of intraocular pressure and glaucoma progression: results from the Early Manifest Glaucoma Trial. Arch Ophthalmol. 2002;120:1268–1279. [CrossRef] [PubMed]
GreaneyMJ, HoffmanDC, Garway-HeathDF, NaklaM, ColemanAL, CaprioliJ. Comparison of optic nerve imaging methods to distinguish normal eyes from those with glaucoma. Invest Ophthalmol Vis Sci. 2002;43:140–145. [PubMed]
TielschJM, KatzJ, QuigleyHA, MillerNR, SommerA. Intraobserver and interobserver agreement in measurement of optic disc characteristics. Ophthalmology. 1988;95:350–356. [CrossRef] [PubMed]
PandaDP, RosenfeldA. Image segmentation by pixel classification in (grey level, edge value) space. IEEE Trans Comp. 1978;27:875–879.
SpenceCD, SajdaP. Role of feature selection in building pattern recognizers for computer-aided diagnosis. Proc SPIE. 2006;3338:1434–1441.
NiemeijerM, StaalJS, van GinnekenB, LoogM, AbràmoffMD. Comparative study of retinal vessel segmentation on a new publicly available database. Proc SPIE. 2004;5370:648–656.
NiemeijerM, vanGinnekenB, StaalJ, Suttorp-SchultenMS, AbràmoffMD. Automatic detection of red lesions in digital color fundus photographs. IEEE Trans Med Imag. 2005;24:584–592. [CrossRef]
LowellJ, HunterA, SteelD, et al. Optic nerve head segmentation. IEEE Trans Med Imag. 2004;23:256–264. [CrossRef]
TailorDR, FinkelLH, BuchsbaumG. Color-opponent receptive fields derived from independent component analysis of natural images. Vision Res. 2000;40:2671–2676. [CrossRef] [PubMed]
SerreT, RiesenhuberM. Realistic modeling of simple and complex cell tuning in the HMAX model, and implications for invariant object recognition in cortex. July 2, 2004;Massachusetts Institute of Technology Cambridge, MA.CBCL Paper 239/AI Memo 2004-004.
QianN, ZhuY. Physiological computation of binocular disparity. Vision Res. 1997;37:1811–1827. [CrossRef] [PubMed]
AbràmoffMD, NiemeijerM. Automatic detection of the optic disc location in retinal images using optic disc location regression. Proceedings of the IEEE EMB Conference on Medical Imaging, New York, NY, August 30–September 3, 2006. 2006;4432–4435.Springer-Verlag New York.
KwonYH, KimYI, PereiraML, MontaguePR, ZimmermanMB, AlwardWL. Rate of optic disc cup progression in treated primary open-angle glaucoma. J Glaucoma. 2003;12:409–416. [CrossRef] [PubMed]
ter Haar RomenyBM. Front End Vision and Multi-scale Image Analysis. 2003;Kluwer Academic Publishers Dordrecht, The Netherlands.
van GinnekenB, Haar RomenyBM. Automatic segmentation of lung fields in chest radiographs. Med Phys. 2000;27:2445–2455. [CrossRef] [PubMed]
KolmogorovV, ZabihR. What energy functions can be minimized via graph cuts?. IEEE Trans Pattern Anal Mach Intell. 2004;26:147–159. [CrossRef] [PubMed]
MarrD. Vision. 2005;Freeman and Company New York.
BirchfieldS, TomasiC. Depth discontinuities by pixel-to-pixel stereo. Proceedings of the Sixth International Conference on Computer Vision, Bombay, India. 1998;1073–1080.IEEE Computer Society Washington, DC.
SonkaM, FitzpatrickJM. Medical Image Processing and Analysis. 2000;The International Society for Optical Engineering Press Wellingham, WA.Handbook of Medical Imaging, vol 2.
SegevR, GoodhouseJ, PuchallaJ, BerryMJ. Recording spikes from a large fraction of the ganglion cells in a retinal patch. Nat Neurosci. 2004;7:1154–1161. [PubMed]
LeeTW, WachtlerT, SejnowskiTJ. Color opponency is an efficient representation of spectral properties in natural scenes. Vision Res. 2002;42:2095–2103. [CrossRef] [PubMed]
OlshausenBA, FieldDJ. Emergence of simple-cell receptive field properties by learning a sparse code for natural images. Nature. 1996;381:607–609. [CrossRef] [PubMed]
AbràmoffMD, KwonYH, Ts’oDY, LiH, BarrigaES, KardonR. A spatial truncation approach to the analysis of optical imaging of the retina in humans and cats. Proc IEEE Int Symp Biomed Imag. 2004;2:1115–1118.
DudaRA, HartPE, StorkDG. Pattern Classification. 2001;Wiley-Interscience, New York New York.
StaalJS, AbràmoffMD, NiemeijerM, ViergeverMA, van GinnekenB. Ridge based vessel segmentation in color images of the retina. IEEE Trans Med Imag. 2004;23:501–509. [CrossRef]
AryaS, MountD, NetanyahuN, SilvermanR, WuA. An optimal algorithm for approximate nearest neighbor searching in fixed dimensions. JACM. 1998;45:891–923. [CrossRef]
ParkinB, ShuttleworthG, CostenM, DavisonC. A comparison of stereoscopic and monoscopic evaluation of optic disc topography using a digital optic disc stereo camera. Br J Ophthalmol. 2001;85:1347–1351. [CrossRef] [PubMed]
CoronaE, MitraS, WilsonM, KrileT, KwonYH, SolizP. Digital stereo image analyzer for generating automated 3-D measures of optic disc deformation in glaucoma. IEEE Trans Med Imaging. 2002;21:1244–1253. [CrossRef] [PubMed]
ScharsteinD, SzeliskiR. A taxonomy and evaluation of dense two-frame stereo correspondence algorithms. ICJV. 2002;47:7–42.
MerickelMB, Jr, WuX, SonkaM, AbràmoffM. Optimal segmentation of the optic nerve head from stereo retinal images. Proc SPIE. 2006;6143:61433B.
CalkinsDJ, TsukamotoY, SterlingP. Microcircuitry and mosaic of a blue-yellow ganglion cell in the primate retina. J Neurosci. 1998;18:3373–3385. [PubMed]
Figure 1.
 
Stereo color image pair of the optic nerve head of the left eye after alignment and cropping to 512 × 512 pixels. The left and right stereo images can be fused by holding the page at approximately 30 cm (12 in.) distance from the eye and gazing into infinity.
Figure 1.
 
Stereo color image pair of the optic nerve head of the left eye after alignment and cropping to 512 × 512 pixels. The left and right stereo images can be fused by holding the page at approximately 30 cm (12 in.) distance from the eye and gazing into infinity.
Figure 2.
 
(AC) Gradings of an optic disc stereo pair by three faculty glaucoma specialists. Rim is grayish, cup is whitish in the left image of the stereo pair from Figure 1 . (D) Reference standard created from (AC) with white, cup; gray, rim; and black, background
Figure 2.
 
(AC) Gradings of an optic disc stereo pair by three faculty glaucoma specialists. Rim is grayish, cup is whitish in the left image of the stereo pair from Figure 1 . (D) Reference standard created from (AC) with white, cup; gray, rim; and black, background
Figure 3.
 
Example of the variety of color opponency, steerable, Gaussian filter bank kernels. Kernels are oriented as if the local gradient were horizontal. First row, from left to right: dark-bright opponency kernels for zeroth order, first order 0° to the local gradient, first-order 90° to the local gradient, second-order 0° to the local gradient, second order 60° to the local gradient, and second order 150° to the local gradient, at scale σ 32 pixels. Second row: same, for scale σ 64 pixels; third row: 128 pixels. Next three rows: same, for blue-yellow opponency kernels; last three rows: same, for red-green kernels. Smaller scales are not shown because they are difficult to see. These kernel images also represent the response of each of the features to an impulse function.
Figure 3.
 
Example of the variety of color opponency, steerable, Gaussian filter bank kernels. Kernels are oriented as if the local gradient were horizontal. First row, from left to right: dark-bright opponency kernels for zeroth order, first order 0° to the local gradient, first-order 90° to the local gradient, second-order 0° to the local gradient, second order 60° to the local gradient, and second order 150° to the local gradient, at scale σ 32 pixels. Second row: same, for scale σ 64 pixels; third row: 128 pixels. Next three rows: same, for blue-yellow opponency kernels; last three rows: same, for red-green kernels. Smaller scales are not shown because they are difficult to see. These kernel images also represent the response of each of the features to an impulse function.
Figure 4.
 
Stereo pair and its most discriminant pixel features. (A, B) The left and right stereo image of the stereo pair. Rows below them show images of the most dominant 12 features in rank order (see Table 1for a description of these features).
Figure 4.
 
Stereo pair and its most discriminant pixel features. (A, B) The left and right stereo image of the stereo pair. Rows below them show images of the most dominant 12 features in rank order (see Table 1for a description of these features).
Figure 5.
 
Classification of stereo pairs by glaucoma faculty, three glaucoma fellows, and pixel feature classification. Rows 1–4: small, medium, and large disc excavation and excavation with inferior notching, respectively. From left to right: left and right images of the stereo pair; reference standard by three glaucoma expert faculty; grading by fellow A, fellow B, and fellow C; and grading by pixel feature classification.
Figure 5.
 
Classification of stereo pairs by glaucoma faculty, three glaucoma fellows, and pixel feature classification. Rows 1–4: small, medium, and large disc excavation and excavation with inferior notching, respectively. From left to right: left and right images of the stereo pair; reference standard by three glaucoma expert faculty; grading by fellow A, fellow B, and fellow C; and grading by pixel feature classification.
Figure 6.
 
Scatterplots of linear cup-to-disc ratio estimates from all 58 stereo pairs by three glaucoma fellows and the pixel feature classification algorithm (vertical axes), against the reference standard grading by three glaucoma faculty (horizontal axes). Clockwise, fellows A and B, the algorithm, and fellow C respectively. r, is the correlation coefficient for each of these estimates with the reference standard.
Figure 6.
 
Scatterplots of linear cup-to-disc ratio estimates from all 58 stereo pairs by three glaucoma fellows and the pixel feature classification algorithm (vertical axes), against the reference standard grading by three glaucoma faculty (horizontal axes). Clockwise, fellows A and B, the algorithm, and fellow C respectively. r, is the correlation coefficient for each of these estimates with the reference standard.
Table 1.
 
Twelve Features Selected as the Most Dominant
Table 1.
 
Twelve Features Selected as the Most Dominant
Rank Description Scale (Pixels) Cumulative Accuracy on Training Set
1 A priori probability of pixel being cup 1 0.810
2 Red-green opponency (opp.) second-order Gaussian derivative transversal to local gradient 32 0.825
3 A priori probability of a pixel being rim 1 0.829
4 Red-green opp. Gaussian blurred (zeroth order) 8 0.832
5 Complex cell output tuned to a disparity of −12 pixels 16 0.845
6 Blue-yellow opp. second-order Gaussian derivative 32 0.854
7 Red-green opp. first-order Gaussian derivative perpendicular to local gradient 16 0.869
8 Bright dark opp. Gaussian blurred (zeroth order) 16 0.872
9 Red-green opp. first-order Gaussian derivative radial to local gradient 4 0.874
10 Blue-yellow Gaussian blurred (first-order) 4 0.874
11 Red-green opp. first-order Gaussian derivative parallel to local gradient 4 0.875
12 Red-green opp. Gaussian blurred (zeroth order) 32 0.876
Copyright 2007 The Association for Research in Vision and Ophthalmology, Inc.
×
×

This PDF is available to Subscribers Only

Sign in or purchase a subscription to access this content. ×

You must be signed into an individual account to use this feature.

×