Automatic cone photoreceptor segmentation using graph theory and dynamic programming

Geometrical analysis of the photoreceptor mosaic can reveal subclinical ocular pathologies. In this paper, we describe a fully automatic algorithm to identify and segment photoreceptors in adaptive optics ophthalmoscope images of the photoreceptor mosaic. This method is an extension of our previously described closed contour segmentation framework based on graph theory and dynamic programming (GTDP). We validated the performance of the proposed algorithm by comparing it to the state-of-the-art technique on a large data set consisting of over 200,000 cones and posted the results online. We found that the GTDP method achieved a higher detection rate, decreasing the cone miss rate by over a factor of five.

To generate quantitative metrics of the photoreceptor mosaic, identification of individual photoreceptors is often a required step. Since manual identification is extremely timeconsuming, many groups have utilized some form of automation when studying the photoreceptor mosaic [9,12,14,17,18,27]. Cone identification algorithms have also been developed and validated for accuracy [39][40][41][42][43]; the Garrioch et al. 2012 algorithm [44], for example, is a modified version of the Li & Roorda 2007 algorithm [39] and was thoroughly validated for repeatability on a large cone mosaic data set. Even so, manual correction was still necessary to identify missed photoreceptors [20,34].
In this work, we propose the use of graph theory and dynamic programming (GTDP), a framework we previously developed to segment layered [45][46][47] and closed contour structures [48], to both identify and segment cone photoreceptors in AO ophthalmoscopy images (Section 2.3). We then validate our algorithm's performance for cone identification (Section 3.2) and evaluate its reproducibility in cone density and spacing estimation (Section 3.3). Finally, the proposed algorithm is extended to segment an image containing both rod and cone photoreceptors (Section 3.4).

Methods
The methods for image acquisition, photoreceptor segmentation, and result validation are discussed in the following sections. Section 2.1 explains the image capture and pre-processing steps, while Section 2.2 describes the gold standard (target) for cone identification. Section 2.3 describes our method for cone segmentation, and Section 2.4 outlines the method for validation. Lastly, Section 2.5 introduces the preliminary rod-cone segmentation algorithm.

Image data set
We validated our algorithm on 840 images (150 × 150 pixels) from the Garrioch et al. study [44], where the methods for image acquisition and pre-processing are described in detail. To summarize, the right eye of 21 subjects (25.9 ± 6.5 years in age, 1 subject with deuteranopia) was imaged using a previously described AOSLO system [13,17] with a 775 nm super luminescent diode and a 0.96 × 0.96° field of view. Four locations 0.65° from the center of fixation (bottom left, bottom right, top left, and top right) were imaged, capturing 150 frames at each site. This process was repeated 10 times for each subject. Axial length measurements were also acquired with an IOL Master (Carl Zeiss Meditec, Dublin, CA) to determine the lateral resolution of the captured images.
Following image acquisition, pre-processing steps were taken in the Garrioch et al. study to generate a single registered image from each 150 image sequence. To do this, first any sinusoidal distortions from the resonant scanner were removed from individual frames. The frames from each sequence were then registered to a reference frame [49], and the top 40 frames with the highest normalized cross correlation to the reference were averaged together. This procedure was performed for all 21 subjects at each of the 4 locations and repeated 10 times over, resulting in a total of 840 images in the image data set. Finally, to ensure that each set of 10 repeated images captured the same patch of retina, the images were aligned using strip registration.
Since the image data set was used strictly for algorithm validation, we obtained a separate set of images to tune the algorithm. These training images were captured using the same imaging protocol, and patients from the test and validation data sets did not overlap.

Gold standard for cone identification
We defined the gold standard as the semi-automatically identified cone locations reported in the Garrioch et al. study, since the cone locations on all 840 images had been carefully reviewed and corrected by an expert grader. As described in the study, the initial cone coordinates were first automatically generated using the Garrioch et al. 2012 algorithm, a modified version of the Li & Roorda 2007 cone identification algorithm [39]. Any missed cones were then added manually. Automatically segmented cones were not removed or adjusted, as the Garrioch et al. 2012 algorithm exhibited a tendency towards false negatives rather than false positives.

GTDP cone segmentation algorithm
We developed a customized implementation of our generalized GTDP framework for closed contour structures [48] to segment cone photoreceptors in AOSLO images. In brief, we used maxima operators to obtain pilot estimates of prominent cones. We then used the quasi-polar transform [48] to map the closed contour cone estimates from the Cartesian domain into layers in the quasi-polar domain. The layered structures were then segmented utilizing our classic GTDP method [45]. By applying the inverse quasi-polar transform, the segmentation lines were carried back into the Cartesian space. Finally, we performed additional iterations to find any missed cones. These steps are described in details in the following.
We first brightened dim photoreceptors by applying Eq.
(1) to the 150 × 150 pixel image indicates a linear normalization of the elements in matrix X to range from y to z.
The range 0.1 to 0.9 was chosen to increase the contrast between the dimmest and brightest pixels, as well as to avoid the log(0) and log (1)   . was determined empirically to avoid thresholding adjacent cones.
To segment each cone, we first used our previously described quasi-polar transform [48] to transform c I to q I (q denotes the quasi-polar domain). To do this, we first transformed c I and c B (Figs. 1(a) and 1(b)) into the polar domain to create p I and p B (Figs. 1(c) and 1(d)).
Next, we column-wise shifted p I until the pilot estimate in p B was flat, resulting in the quasi-polar images q I and q B (Figs. 1(e) and 1(f)). After obtaining , q I we removed regions containing other pilot estimates and already-segmented cones from the search space, and used GTDP to find the shortest path across q I with the following weight scheme: where ab w is the edge weight connecting nodes a and b, LD n g and DL n g are the vertical lightto-dark and dark-to-light gradients [45] of the image at node n, respectively, ab d is the Euclidian distance from node a to b, and . 00001 .
The vertical light-dark gradient comprised the majority of the weight, since it was the primary indicator for the boundary of the central cone. A smaller weight was given to the dark-light gradient to segment boundaries of dimmer cones adjacent to brighter cones ( Fig. 1(c), left). Finally, a vertical distance penalty was added to discourage the segmented line from including adjacent cones. Specific values for weight ranges were determined empirically. We then transformed the shortest path from the quasi-polar domain ( Fig. 1(g)) back into the Cartesian domain to obtain the final segmentation of the cone (Fig. 1(h)), keeping it only if the mean radius was greater than one pixel. This entire process was then repeated for all subsequent cone estimates. At this stage of the algorithm, the cones identified and segmented by the GTDP method ( Fig. 2(b), black) may be similar to those detected by previous methods, since local maxima were used to initialize the cone locations. To further identify any missed cones, we obtained pilot estimates of the cones using a second method: image deblurring using maximum likelihood blind deconvolution [50][51][52] (deconvblind function in MATLAB) with a Gaussian point spread function of half the mean radius of already segmented cones, followed by locating all regional maxima with a pixel connectivity of eight. Any pilot estimates lying outside already-segmented cone locations (Figs. 2(a) and 2(b), white) were segmented using the same quasi-polar GTDP technique, with the modification to the weighting matrix as shown in Eq. (3). In this weighting scheme, the vertical dark-light gradient was assigned a higher weight since cones detected during this section iteration were typically dimmer and adjacent to brighter cones. The vertical distance penalty was also removed since adjacent cones were already segmented and thus removed from the search region.

Statistical validation
We validated our GTDP algorithm by comparing its performance to the Garrioch et al. 2012 algorithm and to the gold standard generated by the Garrioch et al. paper [44]. To perfectly replicate the Garrioch et al. study, all images were cropped to a 55 µm × 55 µm region about the image center to remove any boundary effects.
To evaluate the performance in cone identification, we compared both fully automatic methods (GTDP and Garrioch et al. 2012) to the gold standard using two metrics: # of true positives, those detected by both the fully automatic and gold standard techniques, and # of false positives, those detected by the fully automatic method but not by the gold standard. A cone was considered to be a true positive if it was within a 1.75 µm Euclidian distance from a gold standard cone. This value was chosen since the mean cone spacing reported in the Garrioch et al. study was approximately 3.50 µm; half this value was therefore a reasonable estimate for the cone radius. If an automatically identified cone did not have any gold standard cones within the 1.75 µm distance, then it was tagged as a false positive. Furthermore, more than one automatically identified cone could not be matched to a single gold standard cone, thus yielding the following relationships:  [53]. The reproducibility of each method was assessed by the comparing cone density (number of cones per mm 2 ) and cone spacing (mean distance from each cone to its nearest neighbor) measurements output by each method at each quadrant. The variability in cone density and spacing measurements (characterized by the variance total ) V stemmed from two sources: 1) variability in measurements taken on the same subject, resulting from the method used (within-subject variability; variance ), within V and 2) variability in true values between subjects, resulting from biological variation between subjects (between-subjects variability; variance between ). V Thus, . between within total The reproducibility was characterized using two components: 1) within-subject coefficient of variation (CV), and 2) intra-class (intra-subject) correlation coefficient (ICC). The within-subject CV was defined as the ratio of the square root of within V to the overall mean measurement, where a lower CV indicates a better the method. ICC was defined as the ratio of between V to total , V thus a ratio closer to 1 indicates a better method.

Preliminary GTDP rod-cone segmentation algorithm
To illustrate the potential of this algorithm to segment images containing both rods and cones, we modified the cone segmentation algorithm described in Section 2.3 to segment a rod and cone photoreceptor image (originally 250 × 250 pixels, scaled to 578 × 578 pixels at 0.186 µm/pixel) captured using the new generation of AOSLO systems [17,54]. In this modified version of the algorithm, photoreceptors were segmented with weights determined by Eq. (5), where n i is the intensity of the image at node n, and n r is the distance of node n from the top of the image . Segmentations with radii less than 3.72 µm were considered to isolate rods, and the rest were re-segmented with the weighting scheme in Eq. (6) to isolate cones. The n r distance penalty was removed since cones have larger radii than rods, and the LD n g weights were removed to delineate the prominent hypo-reflective region surrounding cones on AOSLO rather than the high gradient boundary.

Results
Section 3.1 discusses the segmentation results of our method, while Sections 3.2 and 3.3 show quantitative results comparing the performance of our method against the state-of-the-art for cone identification and cone density and spacing reproducibility, respectively. Finally, Section 3.4 shows a preliminary segmentation result for an image containing both rod and cone photoreceptors.

Cone identification performance
The performance in cone identification for each of the methods is shown in Table 1. This table shows that after taking into consideration all correlated data, our GTDP method correctly detected 99.0% of the cones, compared to the Garrioch et al. 2012 method which detected 94.5% of the gold standard cones; this difference was found to be significant (Z = 15.0, p<0.0001). In addition, 1.5% of the cones found by the GTDP method were not in the gold standard. False positive cones could not be detected by the Garrioch et al. 2012 method since the gold standard was based off of the Garrioch et al. 2012 algorithm (see Section 2.2). Lastly, the mean distance error from the true positive GTDP cones to the gold standard cones was 0.20 ± 0.26 µm.  Finally, Fig. 5 takes a closer look at the results from Fig. 4(b) (right). The black box highlights a "false positive" cone added by the GTDP algorithm per the gold standard, however inspection of the original image in Fig. 5(a) indicates that a cone is indeed present at that location. In contrast, the white boxes in Fig. 5 highlight "false negative" cones missed by the algorithm per the gold standard. By inspecting Fig. 5(a), however, these locations do not seem to exhibit hyper reflectivity.  . White boxes: locations where the algorithm "missed" a cone, even though there appears to be no cone present. Black box: location where the algorithm "erroneously added" a cone, although the original image seems to contain an added cone not identified by the gold standard. Table 2 shows the mean, ICC, and within-subject CV values for the cone density and spacing metrics as measured by the Garrioch, GTDP, and gold standard methods separated by image quadrant. The average GTDP cone density ICC of 0.989 indicates that on average, 98.9% of the total variability in the measurements was due to the variability between subjects, while only 1.1% was due to the GTDP algorithm. The average GTDP within-subject CV of 0.0146 indicates that the error in reproducing the same measurement for the same subject was within 1.46% of the mean.

Preliminary rod and cone segmentation result
Figure 6(a) shows an example rod and cone photoreceptor image [17,54] accompanied by the GTDP segmentation result in Fig. 6(b) and its associated centroids in Fig. 6(c). Figure 6(d) shows a histogram of the number of photoreceptors at various sizes based on the segmentation from Fig. 6(b), and Fig. 6(e) demonstrates a simple classification of rod and cone photoreceptors using a size threshold of 27.7 µm 2 .

Discussion and conclusion
We developed a fully automatic algorithm using graph theory and dynamic programming to segment cone photoreceptors in AOSLO images of the retina and validated its performance. We were able to achieve a higher cone detection rate, more accurate cone density and spacing measurements, and comparable reproducibility compared to the Garrioch et al. 2012 algorithm. Furthermore, the segmentation-based approach enabled identification and classification of rods and cones within a single image. This is highly encouraging for largescale ophthalmic studies requiring an efficient and accurate analysis of the photoreceptor mosaic. We obtained the data set from the Garrioch et al. study [44] to validate the performance of our algorithm on a large untrained data set. We compared the performance of our fully automatic cone segmentation algorithm to the state-of-the-art technique, and found that our GTDP method decreased the Garrioch et al. 2012 cone miss rate by a factor of 5.5 (Table 1, 1.0% vs. 5.5% false positives). One point five percent of the cones not identified by the gold standard were also found using our technique. While this implies that our algorithm falsely identified these cones, Fig. 5 shows that in some cases, our GTDP method was able to identify cones not found by the gold standard; such observations, while not the norm, are likely due to the resource intensive nature of semi-automatic cone identification.
The mean results in Table 2 indicate that the cone density and spacing metrics extracted by the GTDP method were more accurate on average than the Garrioch et al. 2012 algorithm, despite the bias where the Garrioch et al. results were used as the starting point for generating the gold standard. While an unbiased comparison could have been conducted, this would have required a fully manual identification of nearly 256,000 cones. Table 2 also shows that the GTDP method generated more reproducible cone density measurements (mean 0.0146 CV) than the other automated method (mean 0.0254 CV). To be consistent with previous publications, we also compared reproducibility in cone spacing, which showed that the Garrioch et al. 2012 method produced more reproducible results (mean 0.0086 CV) compared to both the GTDP method (mean 0.0130 CV) and the gold standard. This is because both the GTDP method as well as the gold standard detected more cones; these were typically the harder and more irregularly spaced cones, and thus resulted in more variable cone spacing. As a result, cone spacing reproducibility might not be the most reliable quantitative measure of performance. Nevertheless, all three methods had a very good within-subject CV, showing that the within-subject standard error (error due to method) ranged by only 0.74% to 2.73% from the mean. Furthermore, all three methods had a very good ICC, showing that 95% to 99.5% of the total variability in the measurements was due to variability between subjects, while only 0.5% to 5% was due to the method. This high ICC was a result of the preprocessing image alignment performed in the Garrioch et al. study (Section 2.1) to ensure that the same patch of retina was imaged.
A notable difference and novelty of the GTDP algorithm as compared to existing en face cone segmentation algorithms, is its use of segmentation to identify cones. While the most common technique for cone identification is to locate points of maximal intensity, such a method only locates cone centers. In contrast, our technique delineates cone boundaries, resulting in added information about the size and shape of the segmented object. This information may be helpful for applications such as studying how the multimodal structure of larger cones changes with time or wavelength. However, it is of importance to note that in the context of AO photoreceptor imaging, cone sizes may be near the resolution limit, especially towards the foveal center. Furthermore, estimation of photoreceptor size depends on the wavelength of the imaging modality (e.g. fundus camera, SLO, OCT) and even varies over time based on intensity fluctuations. As a result, extracting size and shape information about the cones, while helpful, may not be an accurate indication of its true morphologic state.
Another advantage of using segmentation is that it enables a higher cone detection rate. By keeping track of the entire area of a cone rather than only its centroid, we can look for added cones in regions where cones have not yet been found ( Fig. 2(b)). Our technique also provides an advantage for isolating rods and cones within a single image (Fig. 6(e)), as we can readily distinguish between the two types of photoreceptors based on their segmented area in normal retinae. Since accurate photoreceptor classification depends on correctly segmented photoreceptors, however, the rods improperly segmented as cones in Fig. 6(b) resulted in misclassification. A more accurate and robust rod-cone segmentation algorithm moving forward will be essential to improving this preliminary classification result.
A limitation of this study is its rather optimistic validation on higher quality images of normal retina. The AO images taken from diseased retinae, however, are often low in quality and plagued with diverse pathological features. This paper is the first step in introducing a conceptually simple yet robust framework adaptable to incorporating the mathematical and algorithmic innovations necessary for segmenting the more challenging real-world, clinical AOSLO images. Future steps include validation of our rod and cone segmentation algorithm, as well as extension and application of our framework to segment more complicated images of photoreceptors in disease states.