Technical Factors Influencing Cone Packing Density Estimates in Adaptive Optics Flood Illuminated Retinal Images

Purpose To investigate the influence of various technical factors on the variation of cone packing density estimates in adaptive optics flood illuminated retinal images. Methods Adaptive optics images of the photoreceptor mosaic were obtained in fifteen healthy subjects. The cone density and Voronoi diagrams were assessed in sampling windows of 320×320 µm, 160×160 µm and 64×64 µm at 1.5 degree temporal and superior eccentricity from the preferred locus of fixation (PRL). The technical factors that have been analyzed included the sampling window size, the corrected retinal magnification factor (RMFcorr), the conversion from radial to linear distance from the PRL, the displacement between the PRL and foveal center and the manual checking of cone identification algorithm. Bland-Altman analysis was used to assess the agreement between cone density estimated within the different sampling window conditions. Results The cone density declined with decreasing sampling area and data between areas of different size showed low agreement. A high agreement was found between sampling areas of the same size when comparing density calculated with or without using individual RMFcorr. The agreement between cone density measured at radial and linear distances from the PRL and between data referred to the PRL or the foveal center was moderate. The percentage of Voronoi tiles with hexagonal packing arrangement was comparable between sampling areas of different size. The boundary effect, presence of any retinal vessels, and the manual selection of cones missed by the automated identification algorithm were identified as the factors influencing variation of cone packing arrangements in Voronoi diagrams. Conclusions The sampling window size is the main technical factor that influences variation of cone density. Clear identification of each cone in the image and the use of a large buffer zone are necessary to minimize factors influencing variation of Voronoi diagrams of the cone mosaic.


Introduction
Data on the density and packing arrangement of photoreceptors and their normal variance are relevant for evaluating the health of the photoreceptor mosaic and contribute to detecting and characterizing the degradation of mosaic quality due to pathologic alterations occurring in retinal diseases. Several methods have been reported for identification of cones in adaptive optics (AO) retinal images [1][2][3][4][5][6][7][8][9][10][11][12][13]. The precision of each new algorithm was in general assessed by comparing the cone locations with those determined manually or by an algorithm with known accuracy [12][13][14][15][16]. In most cases, careful manual analysis was shown to produce quite accurate cone identification. A moderate to high inter-individual variability (defined as the ratio of the standard deviation to the mean) of cone density, ranging between 12% and 20% between 250 and 1300 mm from the fovea, has been reported in studies of healthy populations [3][4][5][6][7]10,17]. Discrepancies between studies have been related both to demographic and technical factors, such as 1) the inclusion of subjects with different ages or eyes with different axial lengths and refractive corrections, 2) the instrument used for biometry, 3) the model eye used to estimate the retinal image size, 4) the use of individual correction of retinal magnification, 5) the use of foveal center or preferred locus of fixation (PRL) as reference point to define retinal eccentricities and 6) the sampling window area used to count cones. As regards demographic factors, it is well established that A) eyes with longer axial length have lower cone density compared with emmetropic eyes, B) cone density declines with increasing retinal eccentricity, C) density variation is higher close to the fovea than at increasing distance and D) the density of cones along the horizontal meridian is 10% higher than along the vertical meridian at each eccentricity [3][4][5][6][7]10,11,13,14,17]. It is also accepted that there is a tendency towards decreasing density with increasing age [3,14,18].
The majority of studies evaluating the cone mosaic have been carried out using 50650 mm to 64664 mm sampling areas to calculate cone density. This approach was chosen in order to compare data acquired in patients with those shown by Curcio et al. [19][20] in cadaver eyes. On the other hand, recent papers have shown the reliability of using multiple complementary descriptors of the cone mosaic in sampling areas wider than currently used [13,21,22]. This approach has been shown to provide a more comprehensive view of the photoreceptor mosaic geometry. In previous work [15], we evaluated the effect of window size and orientation on cone density at four retinal parafoveal locations, showing that density estimates decrease with decreasing window size. The packing arrangement of cones, estimated using Voronoi diagrams, was not influenced by the sampling window size or its geometry.
The scope of this work was to investigate the influence of technical factors on the variation of cone density and Voronoi diagrams estimated in AO-flood illuminated images of the parafoveal cone mosaic in a population of healthy adults. These technical factors included the sampling window size, the corrected retinal magnification factor (RMF corr ), the conversion from radial to linear distance from the PRL, the displacement between the PRL and foveal center and the manual checking of cone identification algorithm performance.

Methods
All research procedures described in this work adhered to the tenets of the Declaration of Helsinki. The protocol was approved by the local ethical committee (Azienda Sanitaria Locale Roma A, Rome, Italy) and all subjects recruited gave written informed consent after a full explanation of the procedure. Inclusion criteria were an age .18 years old, no history of systemic or ocular diseases and no previous eye surgery. Subjects recruited for the study received a complete eye examination, including non-contact ocular biometry using the IOL Master (Carl Zeiss Meditec Inc, Jena, Germany) and retinal imaging using a Spectralis SLO/SD-OCT (Heidelberg Engineering GmbH, Heidelberg, Germany).
A flood-illuminated AO retinal camera (rtx1, Imagine Eyes, France) was used to acquire images of the cone mosaic. The imaging sessions were conducted after dilating the pupil with one drop of 1% tropicamide. In this study, image sequences of 40 frames each, subtending 4 degrees of visual angle, were recorded at 13 retinal locations, extending between 2 degree nasal, 2.5 degree temporal and 2.5u superior and inferior from the PRL in the right eye of each subject. During imaging, fixation was maintained by instructing the patient to fixate on the internal target of the instrument moved by the investigator.
A proprietary program provided by the manufacturer [23] has been used to correct for distortions within frames of the raw image sequence and to register and frame-average in order to produce a final image with enhanced signal-to-noise ratio. The final AO images acquired at 1.50 degree temporal and superior from the PRL (identified as the point with coordinates x = 0u and y = 0u and here used as the foveal reference point), have been used for subsequent analysis.

Image analysis
We used the nonlinear formula of Drasdo and Fowler and the Gullstrand schematic model eye parameterized by the biometry measurements (corneal central curvature, anterior chamber central depth, axial length) to convert each final image from degrees of visual angle to micrometers on the retina [7,13,15,24,25]. The corrected magnification factor (RMF corr ) was calculated for each eye in order to correct for the differences in optical magnification and thus retinal image size between eyes [13,15,18].
Three areas of different size (3206320 mm, 1606160 mm and 64664 mm) were cropped from each final image and used for subsequent analysis of cone density and preferred packing arrangements of cones. The cone image labelling process was performed using an algorithm implemented with the Image Processing toolbox in Matlab (The Mathworks Inc, Natick MA, USA), as previously described in detail [1,13,15]. Cones were identified independently in each sampling window. The performance of the cone identification algorithm was verified by two expert investigators (ML and GL). When they fully agreed, the position of each cone, that was not automatically labeled, was digitized manually by one investigator by clicking on the cone and marking its location on the image. This procedure ensured that no cone was excluded. Cones whose edges were even partially outside the image section were not labeled. The x, y coordinates of the cones in each sampling window were then stored in a text array and used to calculate the cone density and packing arrangement.  Cone counts were converted into local densities by calculating their number per square millimeter (cones/mm 2 ). Cone density was estimated at both radial and linear distance from the PRL along the horizontal and vertical meridians, which included 1.5 degree and 418 mm temporal and superior from the fovea respectively. The linear distance represented the average distance from the PRL converted from 1.5 degree in our study population. For each image, all the sampling windows were recentered to this linear distance from the PRL. Density values obtained at 1.5 degree eccentricity were also calculated both with and without applying the RMF corr . In the latter case, a single RMF value was used for each sampling area in all subjects. The distance  from the PRL to the foveal center has been further calculated. Although the foveal cones were not resolved by the rtx1 instrument, registration of the montage image of the cone mosaic with the SLO/SD-OCT images helped to identify the foveal center (identified as the foveal pit) in each eye. The cone packing arrangement was analyzed using Voronoi diagrams, as previously described [13,15]. The Voronoi tessellation was implemented by the Matlab voronoi function using the two-dimensional coordinates of the labelled cones. Each Voronoi cell was colour-coded according to the number of neighbouring cones: gray = tetragonal (4n) arrangement, yellow = pentagonal (5n) arrangement, green = hexagonal (6n) arrangement; blue =heptagonal (7n) arrangement and white = octagonal (8n) arrangement. The Voronoi regions containing pixels that extended beyond the bounds of each section were excluded from further analysis, thus creating a thin buffer zone to minimize the boundary effect [26].

Statistics
Retinal data were expressed as mean 6 standard deviation. Statistics were performed using the SPSS software (SPSS Inc., version 17.0). The normal data distribution was verified using the P-P plot within the software. The analysis of variance and the Tukey pairwise test were used to test significance between the cone density measurements and the preferred packing arrangements of cones taken within windows of different size at the same retinal location. The coefficient of variation was used to analyze the variation of cone density between sampling windows of different size at each retinal location.
The intraclass correlation coefficient (ICC; two-way random effects model) was calculated in order to estimate the association between cone density values calculated within the various sampling windows. Bland-Altman analysis [27,28] was used to assess the 95% limits of agreement (LoA) between the cone density values estimated between the various sampling window conditions along the same meridian.
A multiple regression analysis was used to determine the relationships between the eye biometry variables (AxL, SEr) and RMF corr . The statistical significance was set at P,0.05 for all the tests performed.

Results
Fifteen adult subjects (6 males and 9 females) were recruited. The subjects were 21 to 46 years old (29.367.6 years), the spherical equivalent refraction (SEr) ranged between emmetropia and 26.25 D (mean, 22.8762.21 D). The average axial length (AxL) was 24.5661.36 mm (range, 21.66 to 27.04 mm). Normal eye examination was recorded in all cases.
The accuracy of the cone identification algorithm decreased as the sampling window size decreased. The mean percentage of manually identified cones was 0.160.5%, 1.360.7% and 4.963.8% for the 3206320 mm, 1606160 mm and 64664 mm windows respectively.

Influence of corrected magnification factor
The average RMF corr was 0.28260.015 mm/deg (range, 0.258-0.313 mm/deg). The absolute difference in cone density between the sampling windows of the same size with or without the use of individual RMF corr was lower than 120 cones/mm 2 (range, 75-117 cones/mm 2 ; P.0.05) at both retinal locations.

Influence of window size
The average cone density ranged between 3080161313 cones/ mm 2 and 2902362497 cones/mm 2 across the various sampling windows at 1.5 degree temporal location. It ranged between 2823162602 cones/mm 2 and 2612863185 cones/mm 2 at 1.5 degree superior location. The measured cone density decreased with decreasing window size. The absolute difference between cone density estimated in 3206320 mm and 1606160 mm areas were lower than 450 cones/mm 2 at both retinal locations (P = 0.09). The differences were greater (P = 0.07; #2100 cones/ mm 2 ) between density values calculated in 3206320 mm and 64664 mm sampling areas. Table 1 summarizes cone density values and their 95% confidence levels at both retinal locations respectively. Cone density was on average 9% higher along the temporal than the superior location. The coefficient of variation increased from 4.2% to 8.6% and from 9.2% to 12.6% with decreasing window size from 3206320 mm to 64664 mm at the temporal and superior locations respectively.

Influence of retinal conversion distance
The average linear distance from the PRL, corresponding to 1.5 degree eccentricity, was 418619 mm, ranging from 387 to 454 mm. The displacement from 418 mm was statistically signif- icantly correlated to SEr (r = 20.55; P = 0.03) and AxL (r = 0.68; P,0.01): the more myopic and, accordingly, the longer the eye, the higher the RMF corr value. The absolute difference between cone density estimated at radial (1.5 deg) and linear (418 mm) distances from the fovea ranged between 98 and 362 cones/mm 2 . There were no statistically significant differences (P = 0.09) between cone density values calculated within windows of the same size along the same retinal meridian.
The average linear distance from the PRL to the foveal center was 27615 mm, ranging from 11 to 46 mm. The displacement of the PRL from the foveal center was not correlated to SEr (r = 0.18; P = 0.58) or AxL (r = 0.19; P = 0.58).

Correlation and agreement between cone density values
The ICC values between cone density calculated within sampling areas of different size ranged from 0.55 to 0.99 (P, 0.001) and between 0.91 and 0.99 (P,0.001) along the temporal and superior locations respectively (tables 2 and 3). The lowest correlation values were found between cone density calculated in 3206320 mm and 64664 mm windows at the temporal location (ICC,0.70).
A high agreement was found between cone density calculated within the sampling windows of same size with or without using the RMF corr . The agreement between values estimated at radial and linear distance was moderate and the agreement calculated between sampling areas of different size was low. The same results

Regression Model
The multiple linear regression model, including AxL and SEr as predictors, explained 70% of the variance of RMF corr across the population (r = 0.70; P = 0.02), as shown in figure 1. The cone density differences between RMF corr and RMF predicted areas of 3206320 mm and 1606160 mm size were lower than 65-and 40cones/mm 2 at both retinal locations respectively. There were no differences (0 cones/mm 2 ) when comparing data calculated within 64664 mm areas.

Voronoi analysis of the cone mosaic
The percentage of hexagonal Voronoi tiles ranged between 44.563.5% and 47.565.3% at the temporal locations (P = 0.25) and between 49.167.9% and 50.764.6% at the superior locations (P = 0.95). Overall, as the sampling window decreased in size, the standard variation increased. At corresponding retinal locations, the maximum difference (3.1%; P = 0.25) was found between the 5n arrangements calculated between the 3206320 mm and 64664 mm sampling areas at 418 mm temporal from the PRL. The differences between all the other non-hexagonal arrangements were #2% (P.0.05). A summary of the preferred cone packing arrangement calculated at radial and linear distances from the fovea is shown in tables 6 and 7 respectively. Overall, the variation of packing arrangements of cones was influenced by the undersampling effect (i.e., presence of any darks areas across the mosaic), the boundary effect and the manual selection of cones missed by the identification algorithm. Figures 2 and 3 show, in two representative cases, the Voronoi maps created within sampling windows of different size at the temporal and superior retinal locations respectively.

Discussion
The clear identification of cones in a selected area of an AO image of the photoreceptor mosaic represents the key factor to reliably analyze the lattice quality. In addition, it is crucial to develop reliable descriptors to classify the normal distribution and spatial arrangement of the cones. A step that follows is to understand the source(s) of cone density variability in healthy populations of adults; in other words, to ascertain whether cone density variation can be ascribed to demographic or technical/ methodological procedures and which of the above variables can have a greater influence. Knowledge of these factors would permit more accurate evaluation of pathological changes of the cone mosaic. The scope of our work was to understand the weighted influence of various technical factors on packing density estimates of parafoveal cones in AO flood illuminated retinal images of healthy subjects using a commercial device.
Cone density decreased with decreasing window size, as previously shown both in AO-flood illuminated and AO-SLO retinal images of the cone mosaic [15,21]. Density values calculated within 64664 mm areas showed 5-6% and 7-8% lower values than within 1606160 mm and 3206320 mm windows at 1.5 degree temporal and superior locations. The absolute differences between the largest and smallest sampling windows were lower than 2100 cones/mm 2 , approaching statistical significance. The 95% confidence level of cone density calculated in 64664 mm areas did not overlap values obtained across 3206320 mm windows. The intersubject variation in density increased as the window size decreased from 3206320 mm to 64664 mm and density values were 9% lower at the superior than at the temporal Table 6. Preferred packing arrangement of cones (average 6 SD, %) calculated using Voronoi tiles at 1.5 degree temporal and superior eccentricities from the PRL.  retinal location. Cone density estimated in areas of different size showed moderate to low agreement; it is therefore not recommended to compare density between different sampling areas in clinical studies. The error increases when comparing density calculated across 64664 mm (or smaller) sampling areas with values calculated in windows much larger than 64664 mm (e.g., $ 1606160 mm). Although the use of a 3206320 mm sampling area (it is larger than 1 degree) represents an extreme option to analyze the photoreceptor mosaic, it is clear from the present study that the window size should be clarified when comparing data between clinical studies. A moderate agreement was found between density values calculated in sampling areas of the same size at corresponding radial and linear retinal distances. This approach was chosen in order to understand the potential error due to comparing density values calculated at radial or linear distances from the foveal center or PRL in different studies, while assuming an equivalent conversion factor between study populations (e.g., 1 degree = 290 mm). In addition, we found an average distance of 27615 mm between the PRL and the foveal center. The displacement of the PRL to foveal center was not correlated to SEr and AxL. This result was in fair accordance with previous studies using either AO-flood illuminated or AO-SLO instruments [4,29], in which the PRL was found to deviate, on average, 18 mm and 34 mm from the foveal center (identified as the peak cone density) respectively. In previous work [10], we showed that by laterally displacing the center of a 50650 mm sampling window by 18 mm, the displacement error in cone density measurements was 1000 cones/mm 2 and 500 cones/mm 2 at 250-and 1300-mm eccentricity, respectively. This type of error was shown to be lower than the accuracy of the cone identification algorithm. Overall, care should be taken when comparing cone density referred to the PRL or foveal center without considering the proportional error related to the window size (the larger the area, the higher the error) and the retinal eccentricity (the closer to the fovea, the higher the error) in different study populations. Caution should be also taken when comparing cone density measured at radial and linear distance from the foveal reference point without taking into consideration the AxL/SEr of the study populations between different studies.
Both a high correlation (R$0.8) and a high absolute agreement were found between cone density estimated in sampling areas of same size when comparing data calculated with or without using the individual RMF corr . This was the only case for which density data, taken within sampling windows of the same size, could be used interchangeably without incurring errors. Overall, cone density has been demonstrated to be fairly constant in eyes with AxL ranging between 22 and 26 mm in previous studies [3,4,6,7,10,14]. The difference between the RMF corr values calculated using different model eyes has been shown to be small [3][4][5][6][7][8][9][10][11][12][13][14][15]. No direct evidence shows which model eye can be considered the most accurate one [17]. Authors have already used a single RMF corr , ranging from 0.275 to 0.290 mm/deg, to compensate for the retinal distance from 0.6 to 12 degrees eccentricity [3,14,20,24,30]. Drasdo and Fowler showed that the RMF value changes less than 0.02 mm/deg from 0 to 10 degrees at the retina of a schematic eye [24]. In this study, the trend observed for the RMF corr in relation to both AxL and SEr was in agreement with previous reports ( figure 4) [3,4,6,7,10,14,25]. In an effort to translate AO ophthalmoscopy to a clinically valuable tool, it would be desirable to simplify its use, for example avoiding the need to acquire biometry data in each subject/eye and using a single conversion RMF value. According to previous data [1,3,4,6,7,10,14,17,20,24,25] and those from the present work, we can propose to use the following RMF corr values for each SEr range when estimating cone parameters in the parafoveal region The analysis of Voronoi diagrams was shown to provide consistent values between sampling areas of different size; the greatest average difference between 6n tiles, although not statistically significant, was 2.0% and was found between the 3206320 mm and 64664 mm areas at 1.5 degree temporal eccentricity. The average % of non hexagonal arrangements was quite similar between sampling windows of different size (difference #2.4%), except for the 5n tiles across the 3206320 mm and 64664 mm windows at 418 mm eccentricity (average difference of 3.1%; P,0.05). In addition, the % of 6n and 7n tiles was almost stable between sampling areas of different size, while the % of 5n and 8n tended to increase with decreasing window size. It was also valuable to understand that the percentage of 6n Voronoi tiles was higher along the superior (. 4%) than temporal retinal location.
In high-quality AO images of the cone mosaic, abnormalities of the cone mosaic have been shown to occur even when estimates of cone density were within normal limits [21,22,31]. The results from the present work illustrated the reliability of using Voronoi analysis to compare data of 6n preferred packing arrangements taken from sampling areas of different size in the same subject or between subjects [15,32]. Overall, the graphical representation of the cone mosaic geometry is less sensitive to window size than cone density, as previously shown [15]. This is in principle due to the fact that Voronoi tessellation is not dependent on the window size but only on the number of the labelled cones and their relative arrangements. The difference in percentage of the preferred packing arrangements of cones in sampling areas of different size was in part caused by the re-selection of cones in each sampling window. The number of neighbours between some Voronoi tiles of the same cone mosaic changed because of small differences in the position of cones between each image section, as shown in figure 5. Manual re-selection of the unidentified or misidentified cones was revealed to be one of the sources of variation of preferred cone packing arrangements in Voronoi diagrams. The results cannot be explained by imaging artifacts or bias, because we used high-quality images of the same cone mosaic in each subject/eye. However, the number of manually added cones increased with decreasing window size, likely due to relatively increased proportion of dark areas in the image of the cone mosaic (e.g., retinal vessels etc.) [12,13,15]. In order to avoid displacement of point coordinates of the same mosaic and thereby develop reliable Voronoi maps, the cone identification algorithm should segment cone apertures, as recently shown by Chiu et al. [33], and not only their peak intensity. The undersampling and the boundary effects were found to be additional factors that influence Voronoi analysis of the cone mosaic. Presence of dark areas across the image (e.g., any retinal vessels, rods, non-waveguiding cones etc.), leading to ambiguous identification of cones, has been previously shown to decrease the accuracy of Voronoi analysis [15,18,26]. The boundary effect increases as sampling area decreases. This phenomenon implies that a higher proportion of cones are close to the boundary with decreasing window size, so that their nearest neighbours in the effective sampling area may not be their real neighbours in the original population [26]. The use of a large buffer zone, wider than used in the present work, can remove this type of error [15,18,26].
Cone spacing has not been included in the present analysis, because it has been previously shown to be less sensitive to errors in cone labelling than cone density [12,13,34]. Limitations from the present work included the analysis of data obtained from only two retinal locations, although along different retinal meridians. Data obtained at 1.5 degree from the fovea cannot be directly extended to areas closer to the fovea, where cone density is changing most rapidly and therefore differences in cone density estimates between sampling windows of different size could be greater. In addition, the results cannot be generalized to other non commercial AO-flood illuminated systems or to AO-SLO instruments, as discussed previously [13].  2). We revealed the effect of manual re-selection of cones missed by automated identification on variation of the relative arrangements of Voronoi tiles in the same mosaic. Panel A: the central inset shows a 1606160 mm overlapped to a 3206320 mm area. The violet crosses represent the cones that have been identified at exactly the same position in both cases (93%); the red and blue crosses show the cones that have been identified only in 1606160 mm and 3206320 mm area respectively. Panel B: the central inset shows a 64664 mm overlapped to a 1606160 mm area. 72% cones have been identified at exactly the same position. doi:10.1371/journal.pone.0107402.g005 In conclusion, we showed the effect of various methodological factors on cone density and packing arrangement estimates in AO flood-illuminated images of the parafoveal cone mosaic. Cone density estimated within sampling areas of different size cannot be compared in clinical studies. The graphical representation of preferred packing arrangements of cones by Voronoi tiles was slightly affected by window size. Assuming that clear and stable identification of each cone in the image is achieved, Voronoi data obtained from sampling areas of different size (using large buffer zone) at corresponding retinal eccentricity can be compared between different subjects and even in the same subject at different time intervals of observation.