Cone Photoreceptor Cell Segmentation and Diameter Measurement on Adaptive Optics Images Using Circularly Constrained Active Contour Model

Purpose Cone photoreceptor cells can be noninvasively imaged in the living human eye by using nonconfocal adaptive optics scanning ophthalmoscopy split detection. Existing metrics, such as cone density and spacing, are based on simplifying cone photoreceptors to single points. The purposes of this study were to introduce a computer-aided approach for segmentation of cone photoreceptors, to apply this technique to create a normal database of cone diameters, and to demonstrate its use in the context of existing metrics. Methods Cone photoreceptor segmentation is achieved through a circularly constrained active contour model (CCACM). Circular templates and image gradients attract active contours toward cone photoreceptor boundaries. Automated segmentation from in vivo human subject data was compared to ground truth established by manual segmentation. Cone diameters computed from curated data (automated segmentation followed by manual removal of errors) were compared with histology and published data. Results Overall, there was good agreement between automated and manual segmentations and between diameter measurements (n = 5191 cones) and published histologic data across retinal eccentricities ranging from 1.35 to 6.35 mm (temporal). Interestingly, cone diameter was correlated to both cone density and cone spacing (negatively and positively, respectively; P < 0.01 for both). Application of the proposed automated segmentation to images from a patient with late-onset retinal degeneration revealed the presence of enlarged cones above individual reticular pseudodrusen (average 23.0% increase, P < 0.05). Conclusions CCACM can accurately segment cone photoreceptors on split detection images across a range of eccentricities. Metrics derived from this automated segmentation of adaptive optics retinal images can provide new insights into retinal diseases.

PURPOSE. Cone photoreceptor cells can be noninvasively imaged in the living human eye by using nonconfocal adaptive optics scanning ophthalmoscopy split detection. Existing metrics, such as cone density and spacing, are based on simplifying cone photoreceptors to single points. The purposes of this study were to introduce a computer-aided approach for segmentation of cone photoreceptors, to apply this technique to create a normal database of cone diameters, and to demonstrate its use in the context of existing metrics.

METHODS.
Cone photoreceptor segmentation is achieved through a circularly constrained active contour model (CCACM). Circular templates and image gradients attract active contours toward cone photoreceptor boundaries. Automated segmentation from in vivo human subject data was compared to ground truth established by manual segmentation. Cone diameters computed from curated data (automated segmentation followed by manual removal of errors) were compared with histology and published data.

RESULTS.
Overall, there was good agreement between automated and manual segmentations and between diameter measurements (n ¼ 5191 cones) and published histologic data across retinal eccentricities ranging from 1.35 to 6.35 mm (temporal). Interestingly, cone diameter was correlated to both cone density and cone spacing (negatively and positively, respectively; P < 0.01 for both). Application of the proposed automated segmentation to images from a patient with late-onset retinal degeneration revealed the presence of enlarged cones above individual reticular pseudodrusen (average 23.0% increase, P < 0.05).
CONCLUSIONS. CCACM can accurately segment cone photoreceptors on split detection images across a range of eccentricities. Metrics derived from this automated segmentation of adaptive optics retinal images can provide new insights into retinal diseases.
Keywords: nonconfocal split detection, cell segmentation, active contour model, normal database, reticular pseudodrusen N oninvasive imaging of cone photoreceptor mosaics in the living human eye has been enabled by various adaptive optics (AO) ophthalmoscopy modalities. [1][2][3] Quantitative assessment of the mosaic through metrics on AO retinal images, such as cone density and spacing, has shown potential for clinical application 2,4 with substantial efforts already realized toward assembling normative databases. [5][6][7][8][9][10][11][12] To overcome the tedious task of manually identifying cones and to remove the variability of human graders, various automated algorithms have been developed for two types of AO modalities: confocal [13][14][15][16] and nonconfocal [17][18][19] AO scanning light ophthalmoscopy (AOSLO). However, most quantitative metrics have been based on representing each cone as a point. 9 In this work, we focused on using region-based descriptors of cone photoreceptors (i.e., representing cones as a cloud of points as opposed to a single point such as the centroid) as seen by nonconfocal split detection AO, 11 starting with cone diameter. Enlargement of cone photoreceptors has been reported in various diseases. 12,20,21 However, when performed manually, this process is several-fold more time intensive than identification of cone centers. Therefore, we present a novel algorithm that builds upon our previous work 20,22 enabling automated segmentation of cone photoreceptors, and demonstrate its potential value for clinical application.
Cell segmentation is an active area of research in digital pathology and microscopy. 23 Active contour models (ACMs) 24,25 are of particular interest here because of their subpixel accuracy as well as robustness to image noise. [26][27][28][29][30] ACM is a propagation process of deformable and closed contours that is controlled by image forces pulling them to stop at object boundaries. 24,25 If the contour has an explicit mathematical representation, such as a spline, ACM is also called ''snake'' 24 ; if it is implicitly embedded into a high-dimensional function and the actual contour is the zero isosurface of the function, then it is also called ''level set.'' 25 Level sets are often exploited in cell segmentation because they are superior to snakes in handing touching cells. [26][27][28][29] Cell shape priors from a set of manually marked cell contours are often embedded into level-set framework to constrain level-set propagation only near the shape boundary, which can prevent cell oversegmentation due to weak cell boundaries that fail to stop contour propagation. 27,30 However, directly applying these approaches to segment cone photoreceptors is challenging given the anisotropic boundary consisting of dark and bright shading on two opposite sides and little to no contrast on orthogonal sides (Fig. 1A). This three-dimensional appearance with oblique illumination also poses an interesting dilemma on defining the precise location of the boundary, which we will discuss later. To date, there have been very few articles addressing image analysis of split detection images, 17,18,22,31,32 with only one report of cone photoreceptor segmentation 22 (defined as the extraction of the region occupied by cells as opposed to identification of cell centers).
Here, instead of manually creating shape priors, 27,30 we proposed to solve the cone photoreceptor segmentation problem by dynamically establishing circular templates for cone photoreceptors, based on automatic detection of their dark and bright regions, and also cell-specific circular templates, which we call circularly constrained active contour model (CCACM). This improves upon our previous work 22 with a novel snake method to improve contour position estimation by achieving subpixel segmentation. We then demonstrate the potential power of the CCACM approach in streamlining the creation of a normal database of cone photoreceptor diameters across a wide range of eccentricities. Finally, we illustrate the potential clinical utility of the photoreceptor inner segment size as a metric by studying a patient with photoreceptors that appear to be locally enlarged over reticular pseudodrusen.

Data Collection
Research procedures adhered to the tenets of the Declaration of Helsinki and were approved by the Institutional Review Board of the National Institutes of Health. Written informed consent was obtained after the nature of the research and possible consequences of the study were explained to the subjects.
Data were acquired by using a custom-built multimodal AO retinal imager (based on an AOSLO) 11,33 outfitted with a computer-controlled fixation system. 34 Image sequences were corrected for eye motion 35 and manually assembled into montages that included both confocal and split detection images, as previously described. 17 In this study, we imaged 10 subjects with no history of systemic or ocular disease (five female, five male; age range, 22-40 years; mean 6 SD, 26.3 6 5.6 years; additional information in Supplementary Table S1) and also a patient with late-onset retinal degeneration (male, 55 years). For each subject, over 100 retinal locations were imaged, from the fovea out to an eccentricity of approximately 6 mm in the temporal direction. The powers of the 790-and 850-nm light sources measured at the cornea were less than 135 lW and 35 lW, respectively, which were less than the maximum permissible exposure set by the American National Standards Institute standard Z136.1 2014. 36 The retinal scaling factor for conversion from degrees to millimeters was computed by using a paraxial ray trace on a three-surfaced simplified model eye 37 using the subject's biometric information (axial length, corneal curvature, and anterior chamber depth) measured with an IOL Master (Carl Zeiss Meditec, Dublin, CA, USA). White arrows indicate a cone photoreceptor with single extreme region detection; black arrows, cone with more than two region detections. A recovery process finds missing opposing intensity regions to complete the process, allowing the proposed segmentation algorithm to handle cones with single-region detection while still tolerating multiple-region detections. Scale bar: 10 lm.

Automated Cone Segmentation Algorithm
The automated cone segmentation algorithm consists of seven steps ( Fig. 1; detailed mathematical description of the proposed algorithm is provided in the Appendix).
Step 1: Dual Region Detection. Because cone photoreceptors contain both dark and bright regions, multiscale circular voting 17 is used to detect region pairs. Briefly, this algorithm detects gradient magnitude values at the boundaries of circular regions and effectively transfers them into the centers of the local radii of curvature. The region centroid is then determined as the point with the highest response values (Fig. 1B, yellow dots).
Step 2: Dual Region Segmentation. Pairs of dark and bright regions are segmented for the purpose of establishing cell-specific circular templates for segmentation in the subsequent steps. Since dual regions are either dark or bright, their intensity values are either regions of local minima or maxima. Their first image derivatives are therefore close to zero, with second derivatives either positive for dark regions or negative for bright regions. Therefore, we used the Hessian matrix to represent the second image derivatives of the 2D image intensity function (see Appendix).
Dual region detections from Step 1 form seed points to initialize this step. Starting from each seed point, region growing is exploited to include image points with positive or negative Hessian matrix response according to Equation A7 to segment dark and bright regions (Fig. 1C, red and blue regions). Since actual cone segmentation is carried out in a subsequent step, the positioning of this dual region segmentation relative to actual cone boundaries is deemphasized here. Instead, we leverage the robustness of the Hessian matrix response to establish stable circular templates for initializing the actual cone segmentation that follows.
Step 3: Dual Region Connection. This step pairs bright and dark regions corresponding to each cell, according to the following criteria: (1) dark regions are always on the left of bright ones ( Fig. 1C; red and blue regions, respectively), (2) their distance between region detections (Fig. 1B, yellow dots) is less than the expected maximum cone radius (4.5 lm in our datasets, which is consistent with what has been reported in histology 38 ), and (3) if there are multiple region candidates, the one with the largest area is selected (Fig. 1C, black arrow). Dual regions from a single cone photoreceptor can thus be connected (Fig. 1D, green lines). Some cells contained only a single region (Fig. 1C, white arrow). The strategy for recovery of such regions is described in Step 4.
Step 4: Convex Hull Determination. This step determines initial bounding regions for each cell from its dual regions. The morphology of the bounding region is used to determine the size of circular templates. Note that there are usually gaps between dual regions (Fig. 1C). To fill these gaps, two distance functions 39 are computed for each gap starting from the boundary contours of its dual regions, where the distance function measures the shortest distance of each image point to the starting contour. The shortest distance between two contours is used as a filling threshold: the gap is filled with a set of image points if their distance to both contours is less than the shortest contour distance (Equation A10; Fig. 1E, green regions). Finally, the combined dual regions plus filled gaps are imported into the convex hull algorithm 40 to determine the smallest convex polygon (Fig. 1E, black contours).
Recovery Procedure. A recovery procedure is applied to cones that contain only a single region (Figs. 1C-E, white arrows). First, we generate an intensity histogram of all dark regions that were identified from the previous steps. Next, for each cone with a single bright region, we define a trapezoidal arc search area starting from the region detection of the current bright region to find a seed point in the missing dark region. The missing dark region is recovered through a fast matching algorithm 25 to grow the seed point. Following recovery, convex hull algorithm can be used to determine the bounding regions for dual regions, except that there are no filled gaps. In the opposite case (dark regions missing bright regions), the same recovery process is applied except that we invert the intensity values of the original split detection image as I x; y ð Þ ¼ 255 À I x; y ð Þ, where the maximum intensity value is 255. Thus, dark regions are changed to bright, and bright ones are inverted to dark for the purposes of recovery.
Step 5: Circular Template Construction. This step aims to create circular templates from bounding regions to help constrain actual cone segmentation. An ellipse is fit to each bounding region from the previous step, with minor and major axes r 1 and r 2 , respectively. A circular template is constructed centered on the bounding region with radius given by ðr 1 þ r 2 Þ=2, based on the assumption that cone photoreceptors are circular in shape (Fig. 1F, yellow circles). Since these circular templates are constructed dynamically from the image itself, they are already placed in the desired image locations to constrain cone segmentation.
Step 6: Circularly Constrained Active Contour Model. The goal of this step is to find cone boundary contours based on (higher) image gradient magnitudes, which is guided by the circular templates from Step 5. Level-set propagation is used to reduce oversegmentation of any cone photoreceptors that touch (e.g., at lower eccentricities, where they are closer together). The level set is initialized by using the circular template contour, and image forces are defined as the weighted linear combination of image gradient magnitudes and normalized distance from circular template boundaries (Appendix, Equation A16). Level sets are iteratively pulled toward cone boundaries by these image forces (Fig. 1G, cyan contours). Image forces reach a minimum when they approach the boundary balanced by image gradients and the circular template, at which point the iterations are stopped.
Step 7: Active Contour Refinement. Cone boundary contours from CCACM often contain many jagged edges due to pixel space (Fig. 1G, cyan contours), which can discretize the computation of cone photoreceptor diameters. Moreover, cones do not appear to have jagged edges in images (or if they do, it is below the resolution limit of the current system). Therefore, the purpose of this step is to reduce the effects of pixelation and to achieve subpixel segmentation. Instead of using level-set propagation to refine contours, an explicit active contour model (snake) 24 is used because it contains a higher-order contour smoothness term. Enforcing the smoothness term stretches the contour to remove the jagged edges in subpixel space (Fig. 1H, green contours).

Cone Diameter Measurement
The area of each cell contour is calculated by treating it as a polygon and is computed on the basis of this representation. 41 Cone diameter was calculated on the basis of a circle with equivalent area.

Pixel Sampling
Pixel sampling is determined during image acquisition by the field of view (FOV): larger FOVs reduce the number of pixels within each cone, which might cause instability of some image operations in the segmentation algorithm that are computed over pixel space. A commonly used FOV for AO imaging (e.g., 0.46 mm, corresponding to 1.58 of visual angle) can have cones that are only approximately 11 pixels across in size, even at eccentric locations ( Fig. 2). In contrast, smaller FOVs lead to larger numbers of pixels within each cone, but may be less ideal for image acquisition owing to the need to acquire more image locations and the increase in retinal exposure to light when compared to larger FOVs. Therefore, determination of acceptable pixel sampling that results in accurate results is important for cone diameter estimation. Ten retinal regions of interest (ROIs) (60 3 60 lm) from 10 healthy subjects at the eccentricity of approximately 4.5 mm were selected to determine the optimal FOV for cone diameter estimation. Each area was imaged with three different FOVs: 0.23, 0.30, and 0.45 mm. Cone contours were manually marked in the image with FOV of 0.23 mm (Fig. 2, red contours) because it contained the highest pixel sampling. Manually marked contours, including repeated markings from the same expert grader, were compared with the results from CCACM. The segmentation accuracy was evaluated by using the average absolute diameter difference (ADD) and relative diameter difference (RDD) as defined in Table 1. All subsequent results were generated from optimized pixel samplings.

Validation of Cone Segmentation Results
Ten retinal ROIs (60 3 60 lm) from 10 healthy subjects at the eccentricity of approximately 4.5 mm from fixation were selected to evaluate segmentation accuracy. At this eccentricity, approximately half of the image is covered by cones, while the other half is not, providing equal opportunity for false positives and false negatives to appear. These 10 retinal ROIs are different from the ones used for the FOV optimization. All cone photoreceptors from these ROIs were manually segmented by a single expert grader to generate ground truth. For evaluation of segmentation accuracy, in addition to ADD and RDD, we also used average symmetric contour distance (ASD), root mean square symmetric contour distance (RMSD), and maximum symmetric contour distance (MSD) ( Table 1). These metrics are widely used to evaluate accuracy in organ and tumor segmentation. 42 Owing to the ''3D'' appearance of cells, the precise cell boundary can be difficult to define (e.g., whether and how much of the ''shadow'' to include in the cell). To better quantify this issue, for each ROI, five cone photoreceptors were randomly selected (approximately half of the cones that are fully contained within each ROI; Supplementary Fig. S1); we then asked the same grader to repeat manual segmentation on unmarked versions of these cones 8.5 months afterwards (at which point the grader did not retain any memory of the previous markings). This grader did not have any knowledge of the automated segmentation results until the completion of the study. Differences between the first and second contours were analyzed to understand repeatability. We also compared automated segmentation results with the first manual marks, second manual marks, and average of two manual marks.
Let C M and C S be manually marked and automatically segmented contours, containing N 1 and N 2 points, respectively. D M;i and D S;i are ith cone diameters computed from C M and C S . N is the number of cones that is contoured in both C M and C S . D M x; y ð Þ and D S x; y ð Þ are two distance functions computed from C M and C S . ADD is the average absolute difference of cone diameters from manually marked and automatically segmented contours. ASD is the average symmetric contour distance; for each point in the manually marked contour, the closest point in the automatically segmented contour is determined from the Euclidean distance; the same process is performed on each point in the automatically segmented contour; the average of all these shortest distances gives ASD; the value is 0 for a perfectly overlapping contour pair. RDD is the average ratio between absolute diameter difference and diameters from manually marked contours. The MSD metric is also similar to ASD, except that the maximum of all shortest distances is chosen instead of the average. The RMSD metric is similar to ASD; it considers the squared distances between the two sets of contour points, averages the squared values, and then takes the square root.

Quantification of Cone Diameters in Relation to Published Studies
As reported previously, 17 cone photoreceptors cannot currently be reliably resolved by using split detection in the foveal center of healthy retinas. Therefore, at eccentricities ranging from approximately 1.35 to 6.35 mm along the temporal direction, 146 ROIs were selected from 10 healthy subjects (70 3 70 lm ROIs). Whenever needed, the ROIs were slightly shifted to avoid blood vessels. Following automated CCACM segmentation, data were manually corrected to remove erroneous identifications for the purposes of establishing a preliminary normal database. Cone diameters were computed from normal database and plotted against published histology and in vivo data. 11,12,38,43 Cone diameter is a region-based quantitative measurement, which differs from existing point-based measurements such as cone density and spacing. Whether there are relationships between the two remains to be explored. Therefore, we also computed the cone density and spacing for the same 146 ROIs, from a previously developed implementation 9 of the density recovery profile approach. 44,45 Univariate linear regression was performed to understand the correlation between diameter and density as well as diameter and spacing, with 1-way analysis of variance performed to determine statistical significance. Differences between cone diameters measured in males and females were compared by using a 2-tailed, unpooled, paired t-test; axial length and age between males and females were compared by using a 2-tailed, unpooled ttest. Moreover, the coefficient of variation was used to evaluate the relative variability of cone diameters at each ROI by dividing the standard deviation of cone diameters by their mean for each eccentricity.

Demonstration of Value of Cone Diameter Measurements in a Patient With Late-Onset Retinal Degeneration
Cone diameter measurements were performed on images from a patient with late-onset retinal degeneration whose eye contained reticular pseudodrusen lesions 46 that were visible with multimodal imaging, including AOSLO, as has been previously described. 47 Here, we selected four ROIs where single reticular pseudodrusen were clearly visible as hyperreflective regions on confocal reflectance AOSLO across a range of eccentricities (2.85, 3.52, 4.39, and 4.58 mm; each ROI was 175 3 175 lm). CCACM, followed by manual correction, was used to segment cone photoreceptors and measure cone diameters and the coefficient of variation on the selected regions corresponding to lesion (cones directly above individual reticular pseudodrusen) and nonlesion (neighboring cones not directly above individual reticular pseudodrusen). For the manual correction, three expert graders independently identified cones inside ROIs and only cones with agreement between at least 2 of 3 graders were selected (the purpose of this was to ensure agreement of cone identification only, before any segmentation analysis). To calculate cone spacing and density in the nonlesion area, we used ROIs of size 70 3 70 lm. For each lesion, cone diameters in lesion and nonlesion areas were compared by using a 2-tailed, unpooled t-test. Nonlesion cone diameters were compared to the expected normal cone diameter from our data by using a 2-tailed, unpooled, paired t-test. Finally, nonlesion cone spacing and density were compared to expected normal values from histology 48 by using a 2-tailed, unpooled, paired t-test.

Selection of Field of View (Pixel Sampling) for CCACM-Based Cone Segmentation
The average ADD on all 10 ROIs was 0.92, 1.05, and 1.55 lm for the FOVs of 0.23, 0.30, and 0.46 mm, respectively ( Fig. 2; Table 2). The automated segmentation accuracy gradually decreased with increasing FOV, suggesting that the FOV of 0.23 mm is the best choice for the proposed implementation of CCACM. Oversegmentation was the primary source of reduced segmentation accuracy in the case of larger FOVs (Figs. 2B, 2C1), caused by inaccurate dual region segmentation. In addition, the multiscale Hessian matrix requires a sufficient number of image points to be stable (Equation A7). Nevertheless, we found that segmentation accuracy could be near fully recovered by upsampling (bicubically increasing pixel density) the 0.46-mm FOV to 0.23 mm (Fig. 2C2). Interestingly, segmentation accuracy remained stable and sometimes even dropped in the case of further upsampling to a FOV of 0.15 mm or more owing to the presence of inhomogeneous regions inside of dark and bright regions within cones that were emphasized by further upsampling. Upsampling also increased computational cost. In our instrument, the 0.23-mm FOV is the smallest FOV that we use. Therefore, the optimal pixel sampling corresponded to a FOV of 0.23 mm, which was used to perform cone segmentation and diameter estimation in this work.

Accuracy of Cone Segmentation
CCACM accurately segmented cone photoreceptors on 10 subjects with varying levels of image quality: ASD, RMSD, MSD, ADD, and RDD were all within a small fraction of the actual cone size (Fig. 3; Table 3). For the 0.23-mm FOV, the diameter of a cone photoreceptor is approximately 22 pixels. This means that the average contour difference between manually marked and automatically segmented is less than 2 pixels (ASD). The difference in cone diameter is also approximately 2 pixels (ADD).  Manual segmentation of cone photoreceptors was repeatable, with average ASD, RMSD, MSD, ADD, and RDD also within a small fraction of cone size, comparing data from the same grader repeated 8.5 months apart ( Supplementary Fig. S1). These correspond to differences of approximately 1 pixel (ASD) and diameter differences of approximately 2 pixels (ADD). However, for a small subset of cones, there were some notable differences observed (Supplementary Figs. S1D2, S1E2; white arrows), Therefore, we compared our automated segmentation results to both manual segmentations and to the average of two manual segmentations. The difference between CCACM versus average segmentation was close to the difference between two manual segmentations, suggesting that the accuracy of the automated segmentation is close to the achievable accuracy of manual segmentation, but with the benefit of zero variation (since the algorithm will always return the same results for the same data).

Normal Database of Cone Diameters Across Eccentricities
A total of 7441 cone photoreceptors were automatically segmented from 10 subjects across the eccentricities ranging from 1.35 to 6.35 mm along the temporal direction. Erroneously segmented cones (defined as the five types of errors described in the validation section above) were manually excluded to ensure high data quality for the purposes of establishing a normal database (Fig. 4, blue contours). The remaining 5191 cones were selected for inclusion into the normal database. Segmented cones were grouped every 0.3 mm from 1.35 to 6.35 mm (18 bins total). Within each group, the mean and standard deviation of cone diameters and actual eccentricities were computed and plotted against published histology and in vivo data. The resulting automatically segmented cone diameters were in good agreement with existing histology and published data (Fig. 5), except for the results from Andrade da Costa, 43 possibly owing to a speciesdependent variation (Cebus apella) ( Table 4). There is very little human-based data on cone photoreceptor diameters (to our knowledge, the only study is that of Scoles et al. 11 and none that have measured large amounts of cones). Our study was in good agreement with the human histology data from Scoles et al. 11 and was consistent with the two other diameter measurements based on split detection images. 11,12 Interestingly, there was a reduction in cone diameter at the eccentricity of 5.5 to 6.0 mm, consistent with results from Scoles et al. 11 This may potentially be related to an increased presence of rods 48 or also a change in the density of RPE cells 49 at this eccentricity. Although additional subjects are needed to confirm, cone diameters in male subjects were slightly smaller at most eccentricities than those in females ( Supplementary  Fig. S2, P < 0.01), which could not be adequately explained by differences in axial length (male: 24.6 6 1.5 mm, female: 23.   Ten split detection AOSLO images from 10 subjects were manually marked for comparison to automated segmentations. The average computational time of automated segmentation is 1.6 6 0.1 seconds. Additional details are shown in Supplementary Table S2. (Supplementary Table S1). The coefficient of variation was also plotted ( Supplementary Fig. S3, top row).

Relationship Between Cone Size and Packing
There was a statistically significant correlation between cone diameter (averaged across the ROI) with both cone density and spacing (P < 0.001 for both; Fig. 6). This is surprising given the large intersubject variability that is observed in cone density, spacing, and diameter when viewed independently. This intersubject variability is most apparent at eccentricities of 4 to 6 mm ( Supplementary Fig. S4). In the case of cone density versus cone diameter, the inverse quadratic fit was only slightly better than the inverse linear fit (R 2 ¼ 0.45 compared to R 2 ¼   (Table 4). In our results, vertical bars denote the one standard deviation of cone diameters, and horizontal bars standard deviations of eccentricities. † Cone diameters were estimated based on the average of two manually-drawn orthogonal lines.

FIGURE 6.
Correlations between cone diameters with cone density and spacing measured in the same cells. Each dot represents data from one ROI, color-coded for the 10 subjects. Cone diameters represent the average diameters across the ROI. 0.43), which is likely due to the presence of rods at the eccentricities plotted.

Quantification of Cone Swelling in Patient Data
In all four ROIs, cones were significantly enlarged above reticular pseudodrusen lesions, compared to neighboring cones from the same patient (P < 0.05 for each of the four lesions; Fig. 7, Table 5). On average, cones above lesions were 23.0% bigger than their neighboring nonlesion cones. There was no difference between cone diameters measured in nonlesion areas compared to the expected normal cone diameter from our data (P ¼ 0.41). A preliminary assessment of the relative variance of the cone diameters (based on calculation of the coefficient of variation in a sliding window across the image as shown in Supplementary Fig. S3) showed that cone diameters were often more irregular when compared to the expected values from the normal database. The size of the sliding window was matched to that used for the computation of the normal database (70 3 70 lm). We were not able to reliably quantify cone spacing and density for cones above lesions because the individual lesions were too small, and not all cones could be reliably identified owing to image quality issues. Computation of cone density is particularly sensitive to misidentification of even a single cone when the size of the region is small (e.g., lesion areas in Fig. 7D3). This fundamental limitation of cone spacing and density (i.e., that it has to be computed over a contiguous array of cones) illustrates the potential application of cone diameter (which can be computed in a handful of noncontiguous cones). The differences in cone spacing and density between nonlesion cones and expected values from histology 48 were not statistically significant (P ¼ 0.06 and P ¼ 0.16, respectively).

DISCUSSION
We provided the first demonstration of a fully automated region-based cone segmentation algorithm (CCACM) for split detection images of cone photoreceptors. CCACM can quickly and accurately segment cone photoreceptors on split detection AOSLO images through dynamically constructed circular shape priors and active contour propagation. We further demonstrated the utility of using CCACM to efficiently construct a normal database of cone photoreceptor diameters across a wide range of eccentricities, from 1.35 to 6.35 mm, and presented the largest database to date consisting of 5191 cones from 10 subjects. Our results are consistent with previously reported values. Finally, application of this approach to a patient with late-onset retinal degeneration to assess cone photoreceptors above individual reticular pseudodrusen lesions 46 demonstrates the potential clinical utility of this as a potential biomarker for disease. Cone density and spacing are the two most commonly used quantitative metrics for cone photoreceptors, but both suffer from the requirement that a sufficiently large, contiguous patch of cones must be completely and accurately identified. Cone density, in particular, is highly sensitive to misidentification of even a single cone for small ROIs. Although the ROIs above individual reticular pseudodrusen lesions are too small to be used to reliably quantify spacing and density, we illustrated that cone diameter can be measured even on a handful of noncontiguous cones, which is a major advantage for studying data from patients with poor image quality. Here, we showed that even on a lesion-by-lesion basis, we can find patches of cones that are significantly enlarged when compared to their immediate neighbors. In the future, we plan to investigate these photoreceptor-based changes in the context of the underlying RPE (Liu J, et al. IOVS 2017;58:ARVO E-Abstract 304). 50 Automating cone diameter measurements could have important applications for the management, diagnosis, and development and testing of new treatments of many diseases including glaucoma, 51 retinitis pigmentosa, 21 achromatopsia, 52 rod monochromacy, 10 cone-rod dystrophy, 53 and macular telangiectasia type 2 (Scoles DH, et al. IOVS 2014;55:ARVO E-Abstract 5951).
Interestingly, cone diameter was correlated to both cone spacing and cone density. It may be possible to reduce intersubject variability in a dataset of cone spacing or cone density by taking into consideration cone diameters. Generally, one would expect cells to expand to fill the space that is available surrounding them, but studying whether there are natural size constraints to cells relative to their packing could be an important step toward understanding how and why cones become enlarged in disease. In this patient, cones were enlarged by an average of 23.0% (and by as much as 30.5% in one of the lesions). In the future, incorporating rod photoreceptors and Müller cells into a model of size and packing will give a more complete understanding.
Given the ''3D'' appearance of cone photoreceptors in split detection images, we found that it can be challenging to precisely define the location of cell boundaries (e.g., in the presence of ''shadows''). We propose the following guidelines based on our experiences: an object is considered as a cone photoreceptor in a split detection image if (1) it has a ''3D'' appearance with dark and bright regions on both sides and (2) the cone boundary is defined as the locations with the largest intensity gradient magnitude. This second criterion is an important one, since the cone boundaries are not sharp edges.
CCACM is a natural extension of our previous cone identification method, 17 as the framework provides the detections of dark and bright regions to start CCACM. We expect that the integration of CCACM with cone identification algorithms can further improve robustness, particularly for the reduction of false positives. There are limitations of CCACM that can be further improved as well as areas for future work. First, many of the parameters used in CCACM were empirically found to be optimal for a 0.23-mm FOV. Further optimization of parameters throughout the algorithm may lead to additional improvements in performance. Second, incorporation of noncircular templates or other geometric templates may be necessary in cases when cells are significantly elongated. Third, we selected approximately 70% of the automatically segmented cones for inclusion in the normal database, comprising only those cones with the clearest boundaries. Most of the excluded cone photoreceptors were located near the fovea, where the dense packing results in a higher likelihood of ambiguous boundaries (due to boundaries of neighboring cones interacting with each other). Nevertheless, the smaller size of cones and higher density near the fovea meant that a similar number of cones could still be selected for the normal database even with a higher exclusion rate. While we do not expect that exclusion of cones with weak boundaries will affect the accuracy of the normal database, in the future, further enhancement of segmentation accuracy of crowded cones could be achieved by training the cone boundary classifier through deep learning 54,55 and applying the classifier to separate touching cells. 56 Fourth, in addition to cone diameter, second-order metrics based on cone segmentation might lead to additional insights about diseases such as Best disease, in which cones have been reported to be noncircular. 50 Finally, the CCACM approach will also be useful for other recently demonstrated nonconfocal AO images that contain pairs of bright and dark regions, such as ganglion cells imaged by using offset-aperture AO 57 or translucent cells (e.g., horizontal cells) imaged in mice by using knife-edge AO. 58 In conclusion, CCACM provides a novel way to quickly and accurately perform cone segmentation across a wide range of eccentricities. This enabled us to efficiently construct the largest normal database to date of cone photoreceptor diameters, illustrating the potential for CCAM as an analysis tool for quantifying photoreceptor changes that occur owing to disease, which might not be captured by existing approaches.
where p denotes points not at C D or C B , respectively, and w denotes points at C D or C B , respectively. The shortest distance d between C D and C B is calculated as The gap is thus filled (Fig. 1E, where e I is the intensity value that is larger than 80% of the samples in the histogram of dark regions found from previous steps.

Circularly Constrained Active Contour Model (CCACM)
If C q ð Þ : 0; 1 ½ ! R 2 parameterizes a contour, then its Euclidean length is given by The actual cone boundary contour,Ĉ q ð Þ, is set up to stop near (1) circular template boundaries as well as (2) an area that contains large image gradient magnitudes. Image forces are thus composed of these two terms, which achieve a maximum when preparing to stop contour propagation (i.e., the contour velocity should vanish whenĈ q ð Þ is found). The influence of circular templates (Fig. 1F, yellow circles) is introduced by using T x; y ð Þ, which is defined as a normalized distance function (Equation A8) from the template boundary contour C T q ð Þ, given by T x; y ð Þ ¼ D normalized x; y ð Þ; where D normalized x; y ð Þ 2 0; 1 ½ : ðA14Þ If an image point x; y ð Þ is close to C T q ð Þ, T x; y ð Þ is close to 0; otherwise, it is close to 1. Therefore, T x; y ð Þcan be used as the speed term to pullĈ q ð Þ to be close to C T q ð Þ. A second speed term is introduced as M x; y ð Þ, which is motivated by the fact thatĈ q ð Þ should stay in the image points with high image gradient magnitudes.
Note that the ''1'' in the denominator of Equation A15 is a tunable parameter that can be further adjusted to modulate this speed term. Integrating Equations A14 and A15 into Equation A13, the search of cone boundary contourĈ q ð Þ can be mathematically formulated aŝ where j is the curvature of C andÑ is the unit inward normal. Equation A17 shows that the minimization process of Equation 16 is the propagation of active contour C over time step s. To achieve high accuracy and numerical stability, the contour C, as well as its underlying image grids, is typically re-discretized after a few iterations, which is especially critical when the contour C moves a large distance. To alleviate such issue, levelset function / x; y ð Þ is instead used to implicitly describe the active contour C, given by where D x; y ð Þ denotes the Euclidean distance of a point x; y ð Þ to C (Equation A8). Therefore, the contour C can be implicitly represented as C ¼ f x; y ð Þj/ x; y ð Þ ¼ 0g. Such implicit contour representation can avoid the discretization of image grids, and Equation A17 is rewritten as ]/ ]s ¼ F r/ j jdiv r/ r/ j j À rF Á r/; where div Á ð Þ denotes the divergence operator. The initial contour at s ¼ 0 is set to C T q ð Þ, which corresponds to the zero-level set of / at s ¼ 0. Iterative evolution of level-set function / x; y ð Þ over s leads to the identification of cone photoreceptor boundaries.

Active Contour Refinement
Explicit active contour model, also called snake, 24 is given by Here, we set b ¼ 1, f ¼ 1000, and k ¼ 0:2. This emphasizes the smoothness term (especially the second smoothness term describing the contour second derivatives, which constrains the contour according to the thin plate equation to effectively remove the jagged edges). The third image term still plays the role of forcing the contour to stop at the image points with high image gradients. The Euler-Lagrange equation 64  Iteratively evolving C q; s ð Þ over s yields the refined cone boundary contour.

Cone Diameter Estimation
Given a cone boundary contourĈ with N points that is discretely represented as p 0 ; p 1 ; Á Á Á ; p N f g , where p i ¼ x i ; y i ð Þ, i 2 N, and p 0 ¼ p N , the polygon area enclosed byĈ can be computed as 41 Assuming the cone photoreceptor is circular shape, the cone diameter D is calculated as