Automated Photoreceptor Cell Identification on Nonconfocal Adaptive Optics Images Using Multiscale Circular Voting

Abstract

Purpose: Adaptive optics scanning light ophthalmoscopy (AOSLO) has enabled quantification of the photoreceptor mosaic in the living human eye using metrics such as cell density and average spacing. These rely on the identification of individual cells. Here, we demonstrate a novel approach for computer-aided identification of cone photoreceptors on nonconfocal split detection AOSLO images.

Methods: Algorithms for identification of cone photoreceptors were developed, based on multiscale circular voting (MSCV) in combination with a priori knowledge that split detection images resemble Nomarski differential interference contrast images, in which dark and bright regions are present on the two sides of each cell. The proposed algorithm locates dark and bright region pairs, iteratively refining the identification across multiple scales. Identification accuracy was assessed in data from 10 subjects by comparing automated identifications with manual labeling, followed by computation of density and spacing metrics for comparison to histology and published data.

Results: There was good agreement between manual and automated cone identifications with overall recall, precision, and F1 score of 92.9%, 90.8%, and 91.8%, respectively. On average, computed density and spacing values using automated identification were within 10.7% and 11.2% of the expected histology values across eccentricities ranging from 0.5 to 6.2 mm. There was no statistically significant difference between MSCV-based and histology-based density measurements (P = 0.96, Kolmogorov-Smirnov 2-sample test).

Conclusions: MSCV can accurately detect cone photoreceptors on split detection images across a range of eccentricities, enabling quick, objective estimation of photoreceptor mosaic metrics, which will be important for future clinical trials utilizing adaptive optics.

Adaptive optics scanning light ophthalmoscopy (AOSLO),^1–6 a type of adaptive optics (AO) ophthalmoscopy,^7–9 has enabled the visualization of individual retinal cells, such as cone photoreceptors, directly in the living human eye. Analyzing the state of the cone photoreceptor mosaic during disease onset and progression is important not only for understanding the cellular nature of retinal diseases, but also for more rapidly evaluating the efficacy of treatments.^9,10 Thus far, the majority of cone photoreceptor imaging has been performed using confocal AOSLO^1–6 or flood-illumination AO^7,11,12 in which single scattered light is captured. To date, relatively few image analysis tools have been developed for AOSLO modalities that utilize multiple scattered light,^13–15 such as nonconfocal split detection.¹⁶ In this paper, we will focus specifically on nonconfocal split detection,¹⁶ which is distinct from split detection that includes the confocal portion of the signal.¹⁷ Split detection is particularly valuable for cone analysis in eccentric locations because the corresponding confocal reflectance images (hereinafter referred to as “confocal”) often contain areas where the presence of cones is ambiguous (Fig. 1). This ambiguity arises from the possibility of not always having a one-to-one correspondence between reflections from cones in the confocal channel and their corresponding cones visible in the split detection channel (white circles, Figs. 1D, 1H), or from the presence of rod photoreceptors, which could be misidentified as cone photoreceptors (white circles, Figs. 1B, 1F). It is known that cone photoreceptors sometimes appear dark on the confocal modality,^18–21 which can hinder their identification (white circles, Fig. 1C, 1G). Thus, split detection images are typically preferred over confocal images for eccentric cone identification. The limitation of current split detection imaging is that cone photoreceptors are difficult to distinguish at or near the fovea (Fig. 1E). This work aims to develop an approach for computer-aided analysis of cone photoreceptors on split detection AOSLO images. 

Figure 1

 

Confocal reflection (top) and split detection (bottom) AO images captured simultaneously at eccentricities of 0.06, 1.5, 3.0, and 4.5 mm from subject 7. Cone density decreases with increasing eccentricities (left to right). Cone photoreceptors are more easily identifiable on the confocal image near the fovea in (A) when compared to its corresponding split detection image (E). In contrast, cone photoreceptors are more reliably identifiable at eccentric locations in split detection images (white circles). Challenges for automated cone identification on split-detection images include irregular light illuminance (white arrow, F), opposite intensity extremes within a cone (white arrow, G), and vascular shadows (white arrow, H). Scale bar: 50 μm.

A prerequisite for quantifying cone photoreceptor mosaic metrics such as density and spacing is accurate cone identification, which is both tedious and time-consuming when performed manually. Existing approaches to automate this step have been mainly focused on image analysis on confocal and flood-illumination AO images, based on local maximum intensity value identification,²² graph theory, and dynamic programming,²³ circle Hough transform,²⁴ or automated detection of the radius of Yellott's ring.²⁵ Unfortunately, differences in the appearance of split detection and confocal modalities prevent the direct application of these existing approaches to nonconfocal images. Individual cones tend to have an intensity pattern resembling a Gaussian profile in confocal AO images (Figs. 1A–D). However, in the case of split detection, the intensity pattern corresponding to a single cell is more similar to Nomarski differential interference contrast microscopy²⁶ in which pairs of dark and bright regions are present in a consistent orientation across the entire image (Figs. 1E–H). In this paper, we demonstrate an automated algorithm to identify cone photoreceptors on split detection images. This identification algorithm is the first step toward automated computation of additional structural descriptive metrics²⁷ for assessing the cone mosaic. 

Cone photoreceptors visualized using split detection AOSLO images have several unique features that can hinder automated cone identification, including irregular light illuminance of the local surrounding tissue (arrow, Fig. 1F), opposing intensity extremes within a cone (arrow, Fig. 1G) and anisotropic contrast at cell boundaries. Moreover, image regions at vascular shadows (arrow, Fig. 1H) often contain unreliable image information. Finally, variations in the size and density of individual cones depending on the eccentricity can require an adaptive algorithm to handle these differences, which is known as the scale-selection problem in the field of computer vision.²⁸ In our previous work utilizing AO images from multiple eccentricities, we showed that the optimal scale can be determined locally according to multiscale circular voting (MSCV).²⁹ In this paper, we further extend upon this MSCV approach²⁹ by: (1) introducing a new MSCV algorithm for identifying and connecting dark and bright region pairs within each cone photoreceptor; (2) refining cone positions by incorporating size estimations derived from a multiscale strategy for region size determination; and (3) extending the identification results to quantitative analysis of density and spacing across different retinal eccentricities with comparisons to histology²⁶ and existing published data.^27,30–33 

Methods

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. 

A custom-assembled AOSLO with split detection capability^2,16 was used to image 10 subjects with no history of systemic or ocular disease (three female, seven male; age range, 20–41 years; mean ± SD, 28.3 ± 6.8 years; additional information in Supplementary Table S1). Two wavelengths were used: 790 nm for retinal imaging and 850 nm for wavefront sensing. The powers measured at the cornea were less than 130 and 25 μW, respectively, which were less than the maximum permissible exposure set by the American National Standards Institute standard Z136.1 2014.³⁴ Prior to imaging, both eyes were dilated with 2.5% phenylephrine hydrochloride and 1% tropicamide. The eye that subjectively yielded better image sharpness in the confocal imaging channel was selected at the start of each imaging session. A computer-controlled fixation system³⁵ was used to assist with the video acquisition at various overlapping retinal locations from the fovea out to an eccentricity of approximately 6 mm in the temporal direction. Videos were acquired using multiple square fields of view (FOV) corresponding to approximately 0.23, 0.30, and 0.45 mm on the retina. The largest FOV was used to more efficiently connect all imaging areas with sufficiently overlapping areas, important for montaging the collected videos for determination of retinal eccentricity. Smaller FOVs, which were acquired at selected regions and registered to the larger FOVs, were utilized for cone identification because they contained higher sampling. We recruited 5 out of 10 subjects for a short acquisition protocol (less than 50 videos total, acquired near the fovea and at select eccentric retinal locations); the remaining five were recruited for an extended acquisition protocol (over 100 videos, acquired near the fovea and along a continuous strip of retinal eccentricities). In order to capture a more representative dataset containing varying levels of image quality, the first 10 subjects to successfully complete either the short or long imaging protocols were used for this study. Since the primary goal of this study was to demonstrate the capabilities of the MSCV algorithm, a sample size of 10 subjects was deemed sufficient and the power of this study was not calculated (to reach the desired sample size, 12 subjects were recruited, with two excluded due to unsuccessful imaging sessions: one due to excessive eye motion and the other due to an inability to sit still due to frequent coughing). 

Raw videos were preprocessed by correcting for eye motion,³⁶ and then averaged and assembled into montages that included both confocal and split detection images. For the purposes of measuring retinal eccentricity, the fovea was manually determined based on a subjective search for the retinal region with the highest cone density on the confocal images. The retinal scaling factor for conversion from degrees to millimeters was computed using a paraxial ray trace on a three surfaced simplified model eye³⁷ using the subject's biometric information (axial length, corneal curvature, and anterior chamber depth) measured using a commercial device (IOLMaster; Carl Zeiss Meditec, Dublin, CA, USA). 

Automated Cone Identification

An automated cone identification algorithm was developed, consisting of four main steps: multiscale circular voting, region size determination, point pair connection, and centering. Detailed descriptions of the algorithm are provided in the Appendix. 

Step 1: Multiscale Circular Voting

The goal of this step is to detect bright and dark regions within cone photoreceptors (blue and yellow crosses in Fig. 2B) by identifying circular objects using a combination of iterative circular voting³⁸ and a multiscale strategy which permits variations in the sizes of cone photoreceptors (Fig. 1). 

Figure 2

 

Overview of cone photoreceptor identification algorithm on split detection AOSLO images. (A) Original image; (B) MSCV of dual intensity extreme regions in each individual cell (blue and yellow crosses); (C) point pair connections (orange lines); (D) region size determination, denoted by red circles; (E) identification centering; (F) final cone identifications (green crosses). The white arrow indicates a cone photoreceptor without a point pair connection. Scale bar: 10 μm.

MSCV begins with an iterative, coarse-to-fine search process. First, for a given starting point with large image gradient magnitude (Fig. 3A), an angular search range of the voting area is cast (white lines in Figs. 3B–E). The actual voting area is a trapezoidal arc (Equation A.3, Fig. A1). The MSCV response is computed within the trapezoidal voting area by convolving image gradient magnitudes with a Gaussian function (Equation A4). Wherever there is radial symmetry (as in the case of circular cells), MSCV responses from different voting areas accumulate in the centers of concave objects (gray/white spots in Figure 3B–E; see also Fig. A2). Note that the MSCV responses correspond to centers of bright and dark regions within photoreceptor cells, which are different than the actual center of the photoreceptor cell, which will be addressed in Step 4. By iteratively shrinking the angular search range, the MSCV response gradually converges to the region centers (Fig. 3E). 

Figure 3

 

Process of MSCV. From left to right: original image, followed by MSCV response image of bright regions after each iteration (top row); their corresponding cone identifications (bottom row). Note that MSCV response gradually converges to the cell region center due to coarse-to-fine strategy controlled by decreasing the voting angle range at cell boundaries after each iteration. Meanwhile, the number of false positives gradually drop, and true identifications are eventually found. The black dot in (A) represents an example boundary point at which the algorithm begins, with white lines denoting the shrinking voting ranges at various steps within the iteration. White x's are image points with maximal MSCV responses in each iteration pertaining to the voting area that is marked with the white lines. Scale bar: 10 μm.

Next, this entire process is repeated across different scales (Equation A1). Here, it is important to determine what a reasonable radius is for the angular search range. If it is too small, the MSCV response at the region center cannot accumulate enough high image gradient magnitudes; if it is too large, image information from other cells will be improperly attributed to each cell's response calculations. To search for a reasonable radius, we determine the MSCV response across a sequence of search radii, and then search for the radii that result in the maximal MSCV response. Specifically, the sequence of search radii generates point clusters which are then converted into a connected image patch using the morphologic dilation shape operator.³⁹ Maximal MSCV responses are then determined for each connected patch (generally each patch corresponds to a single cone photoreceptor). The radius that generated this maximal MSCV response is then selected and assumed to be reasonable. 

Step 2: Point Pair Connection

This step connects the bright and dark regions within each of the cells, according to the following criteria: (1) dark region centers are always on the left of bright ones (yellow and blue crosses, respectively in Fig. 2C); and (2) their distance is less than the expected maximum cone radius (4.5 μm in our datasets which is consistent with what has been reported in histology⁴⁰). If there are multiple possibilities for region pairs, the one with the smallest distance is selected. Some cells contained only a single detection (white arrow, Fig. 2C; approximately 20% of all cells in our validation datasets). The strategy for the identification of unpaired cells is described in Step 4. 

Step 3: Region Size Determination

In addition to identifying the dual region centers, the multiscale strategy is also important for determining the sizes of each region, which are used in the subsequent step to determine the final positioning of each identification. This step utilizes the local image structure and a scale-invariant image operator to perform optimal scale selection (based on scale-space theory).²⁸ 

Scale selection is initiated by constructing a one-parameter family of smoothed images based on a scale parameter (Equation A8). Within this multiscale representation of the original image, the scale parameter is treated as the additional continuous variable. According to the scale selection theory,²⁸ the optimal scale is determined by searching for local extrema of a scale-invariant image operator such as the scale-normalized Laplacian operator^41,42 (Equation A10). This step outputs circles representing the scale of each identified region which is approximately equivalent to its size (Fig. 2D). 

Step 4: Centering

The final position of each identified cone photoreceptor is determined by taking the weighted average of their point pairs from Step 2, where weights are defined as their optimal scales from Step 3. By using a weighted average to determine cell centers (Equation A11) the final identification can be placed closer to the centers of each cell (Fig. 2F). For any unpaired cells (white arrows, Fig. 2), we only keep single identifications from bright regions. In this case, the scale-normalized Laplacian operator is used to locally adjust the locations of the single identifications by moving positions until Laplacian value achieves the local extrema in the multiscale image space, as demonstrated in our previous approach.²⁹ 

Quantitative Analysis of Cone Photoreceptor Cell Identifications

Averaged split detection images were assembled into a montage of overlapping retinal regions. The montage was used only to guide the selection of regions of interest (ROI) for cone mosaic analysis at different eccentricities. In order to avoid distortion artifacts due to the use of overlapping images that were calculated from different reference frames,⁴³ each ROI was fully contained within a single image of the montage. The sizes of ROIs were set to 55 × 55 μm for eccentricities less than 3 mm and 100 × 100 μm otherwise. Voronoi diagrams were constructed for each ROI using detected points except for those that generated Voronoi neighborhoods that exceeded the boundary of the ROI. Cone density was then computed by dividing the number of cones by the total area of the Voronoi neighborhoods.⁴⁴ Cone spacing was estimated according to a previously developed implementation²⁷ based on the density recovery profile.^45,46 

Validation of Identification Accuracy Across Different Eccentricities

Thirty images were selected at the temporal eccentricities of 1.5, 3.0, and 4.5 mm from all 10 subjects (Supplementary Table S1). These eccentricities were selected to evaluate the robustness of MSCV, since cone density varies substantially across these three eccentricities (Fig. 1). For validation purposes, ROIs of 440 × 440 pixels were generated for all images. To generate ground truth, cone photoreceptors were manually labeled by the same grader in all selected images. In order to evaluate the accuracy of automated identifications in comparison to ground truth, recall, precision, and F1 score were computed (commonly used for the evaluation of cell identification accuracy⁴⁷). These metrics range from 0 to 1, where low values indicate poor matching between automatic cone identifications compared to ground truth, and high values indicate good matching.  

\(\def\upalpha{\unicode[Times]{x3B1}}\)\(\def\upbeta{\unicode[Times]{x3B2}}\)\(\def\upgamma{\unicode[Times]{x3B3}}\)\(\def\updelta{\unicode[Times]{x3B4}}\)\(\def\upvarepsilon{\unicode[Times]{x3B5}}\)\(\def\upzeta{\unicode[Times]{x3B6}}\)\(\def\upeta{\unicode[Times]{x3B7}}\)\(\def\uptheta{\unicode[Times]{x3B8}}\)\(\def\upiota{\unicode[Times]{x3B9}}\)\(\def\upkappa{\unicode[Times]{x3BA}}\)\(\def\uplambda{\unicode[Times]{x3BB}}\)\(\def\upmu{\unicode[Times]{x3BC}}\)\(\def\upnu{\unicode[Times]{x3BD}}\)\(\def\upxi{\unicode[Times]{x3BE}}\)\(\def\upomicron{\unicode[Times]{x3BF}}\)\(\def\uppi{\unicode[Times]{x3C0}}\)\(\def\uprho{\unicode[Times]{x3C1}}\)\(\def\upsigma{\unicode[Times]{x3C3}}\)\(\def\uptau{\unicode[Times]{x3C4}}\)\(\def\upupsilon{\unicode[Times]{x3C5}}\)\(\def\upphi{\unicode[Times]{x3C6}}\)\(\def\upchi{\unicode[Times]{x3C7}}\)\(\def\uppsy{\unicode[Times]{x3C8}}\)\(\def\upomega{\unicode[Times]{x3C9}}\)\(\def\bialpha{\boldsymbol{\alpha}}\)\(\def\bibeta{\boldsymbol{\beta}}\)\(\def\bigamma{\boldsymbol{\gamma}}\)\(\def\bidelta{\boldsymbol{\delta}}\)\(\def\bivarepsilon{\boldsymbol{\varepsilon}}\)\(\def\bizeta{\boldsymbol{\zeta}}\)\(\def\bieta{\boldsymbol{\eta}}\)\(\def\bitheta{\boldsymbol{\theta}}\)\(\def\biiota{\boldsymbol{\iota}}\)\(\def\bikappa{\boldsymbol{\kappa}}\)\(\def\bilambda{\boldsymbol{\lambda}}\)\(\def\bimu{\boldsymbol{\mu}}\)\(\def\binu{\boldsymbol{\nu}}\)\(\def\bixi{\boldsymbol{\xi}}\)\(\def\biomicron{\boldsymbol{\micron}}\)\(\def\bipi{\boldsymbol{\pi}}\)\(\def\birho{\boldsymbol{\rho}}\)\(\def\bisigma{\boldsymbol{\sigma}}\)\(\def\bitau{\boldsymbol{\tau}}\)\(\def\biupsilon{\boldsymbol{\upsilon}}\)\(\def\biphi{\boldsymbol{\phi}}\)\(\def\bichi{\boldsymbol{\chi}}\)\(\def\bipsy{\boldsymbol{\psy}}\)\(\def\biomega{\boldsymbol{\omega}}\)\(\def\bupalpha{\bf{\alpha}}\)\(\def\bupbeta{\bf{\beta}}\)\(\def\bupgamma{\bf{\gamma}}\)\(\def\bupdelta{\bf{\delta}}\)\(\def\bupvarepsilon{\bf{\varepsilon}}\)\(\def\bupzeta{\bf{\zeta}}\)\(\def\bupeta{\bf{\eta}}\)\(\def\buptheta{\bf{\theta}}\)\(\def\bupiota{\bf{\iota}}\)\(\def\bupkappa{\bf{\kappa}}\)\(\def\buplambda{\bf{\lambda}}\)\(\def\bupmu{\bf{\mu}}\)\(\def\bupnu{\bf{\nu}}\)\(\def\bupxi{\bf{\xi}}\)\(\def\bupomicron{\bf{\micron}}\)\(\def\buppi{\bf{\pi}}\)\(\def\buprho{\bf{\rho}}\)\(\def\bupsigma{\bf{\sigma}}\)\(\def\buptau{\bf{\tau}}\)\(\def\bupupsilon{\bf{\upsilon}}\)\(\def\bupphi{\bf{\phi}}\)\(\def\bupchi{\bf{\chi}}\)\(\def\buppsy{\bf{\psy}}\)\(\def\bupomega{\bf{\omega}}\)\(\def\bGamma{\bf{\Gamma}}\)\(\def\bDelta{\bf{\Delta}}\)\(\def\bTheta{\bf{\Theta}}\)\(\def\bLambda{\bf{\Lambda}}\)\(\def\bXi{\bf{\Xi}}\)\(\def\bPi{\bf{\Pi}}\)\(\def\bSigma{\bf{\Sigma}}\)\(\def\bPhi{\bf{\Phi}}\)\(\def\bPsi{\bf{\Psi}}\)\(\def\bOmega{\bf{\Omega}}\)\begin{equation}\tag{1}Recall = {{TP} \over {TP + FN}},Precision = {{TP} \over {TP + FP}},{\rm{F}}1 = {{2TP} \over {2TP + FP + FN}}\end{equation}

Here, TP refers to the true positive of cone identification; FP, false positive; and FN, false negative.  

Image quality plays an important role on the identification accuracy. However, there are currently no strategies for objectively evaluating the image quality of split detection images. Here, we propose the use of anisotropy to objectively evaluate image quality.⁴⁸ In brief, this method assumes that high-quality images contain dominant directions of information content, and therefore applies a pseudo-Wigner distribution to measure a pixel-wise entropy along different directions. Anisotropy is defined as the standard deviation of entropy values (normalized from 0 to 1). Larger anisotropy values indicate the presence of dominant directions interpreted to mean higher-quality images. We used this approach to objectively assess split detection images for interpretation of our accuracy results, based on the assumption that high-quality split detection images have high contrast along the horizontal direction. 

Validation of Spacing and Density Metrics Computed From Automated Identifications

Cone photoreceptors were difficult to visualize near the fovea in split detection images (Fig. 1E). Thus, 146 ROIs were selected from 5 subjects undergoing the extended-period protocol at eccentricities ranging from approximately 0.4 to 6.0 mm along the temporal direction (Supplementary Table S1). ROIs were used to compute cone density and spacing based on MSCV-derived cone identifications for direct comparison to previously published adaptive-optics-based measurements in normal subjects (Table 1). In addition, the relative error between MSCV-derived measurements and histology was computed by first averaging both density and spacing measurements every 0.3 mm and then by computing the relative error, Display Formula\(\varepsilon \).  

\begin{equation}\tag{2}\varepsilon = \mathop \sum \limits_i^N {{\left| {v_i^a - v_i^h} \right|} \over {v_i^h}}\end{equation}

Here, Display Formula\(i\) refers to the index of selected eccentricities, Display Formula\(v_i^a\) is the averaged MSCV-based cone density or spacing at eccentricity corresponding to the index Display Formula\(i\), and Display Formula\(v_i^h\) is the corresponding histology value. The Kolmogorov-Smirnov 2-sample test was also used to evaluate whether there were any differences between MSCV-based density measurements and histology.  

Table 1

 

AO-Based Normal Databases for Cone Density and Spacing

 

Comparison Between Automated Cone Identification From Confocal and Split Detection Images

Five pairs of simultaneously acquired confocal and split detection images were selected from five subjects and independently analyzed using automated identification algorithms (Supplementary Table S1). For each subject, the largest ROI that contained an unambiguous cone mosaic and relatively constant cone density on both the confocal and the split detection channels was selected, which was empirically 400 × 400 pixels in size at an eccentricity of about 0.75 mm (average ± standard deviation: 0.73 ± 0.09 mm). At this location, confocal images contain relatively few rods and split detection images start to have clear cell boundaries (between the eccentricities shown in Figs. 1A, 1E and Figs. 1B, 1F). A previously developed automated cone identification algorithm was used to analyze the confocal images²² for direct comparison to MSCV-based analysis of the corresponding split detection image. Recall, precision, and F1 score were again used to evaluate consistency. We also compared cone spacing and density measurements from two image modalities. 

Results

High Identification Accuracy Across Different Eccentricities

MSCV accurately detected cone photoreceptors in five subjects with average recall, precision, and F1 score of 92.9%, 90.8%, and 91.8% (Fig. 4; Table 2). The average computational time was 0.9 seconds on the image of 440 × 440 pixels (Supplementary Table S2; overall, it was faster with increasing eccentricities). On average, there was no clear difference in accuracy across the three eccentricities tested. Identification accuracy was higher in the case of better image quality (subject 2, Figs. 4A–C) but lower in the subject with the worst image quality, due to a higher rate of false negatives (subject 3, Figs. 4D–F). Image quality assessment based on anisotropy also confirmed that detection accuracy was more likely to be higher on images with larger assessment values (Supplementary Figs. 1, 2). False negatives were likely due to weak cell boundaries (yellow marks in Figs. 4E, 4F). False positives were mostly attributed to either areas of unusually high brightness leading to multiple identifications within one region (orange circle in Fig. 4C), and/or irregularities in the areas between cones, possibly caused by a weak rod photoreceptor signal (orange circle in Fig. 4E). Artificial boundaries created by multiple nearby photoreceptors in close proximity to each other also gave rise to both false negatives and positives (yellow marks in Figs. 4A, 4D). 

Figure 4

 

Identification results comparing automated MSCV identifications to manual labeling for subject 2 (A–C) and subject 3 (D–F) at eccentricities of 1.5 (A, D); 3.0 (B, E); and 4.5 mm (C, F). True identifications, green crosses; false positives, orange crosses; false negatives, yellow triangles. Examples of some typical false positives (multiple bright region detections within a single cell) and false negatives (cone photoreceptors with weak boundaries that are not identified) are highlighted with circles. Scale bar: 50 μm.

Table 2

 

Evaluation of Identification Accuracy and Computational Time

 

Reliable Estimation of Quantitative Cone Photoreceptor Mosaic Metrics

Results were in good agreement with both histology²⁶ and published data^27,30–33 in the eccentricity range tested here (0.4–6.0 mm, Fig. 5). There was no statistically significant difference between the histology and MSCV-based density measurements (P = 0.96, Kolmogorov-Smirnov 2-sample test⁴⁹). The relative errors between MSCV-based cone density and spacing measurements in comparison to histology were 10.7% ± 6.4% and 11.2% ± 4.8%, respectively. These relative errors can be largely attributed to individual-to-individual variations, which would likely average out with larger sample numbers (Fig. 6). Exponential fitting was performed on the MSCV-based cone spacing results as has been previously described according to Equation 3:^50,51  

\begin{equation}\tag{3}spacing = A \times {e^{ - B \times eccentricity}} + C\end{equation}

where Display Formula\(A,B\) and Display Formula\(C\) are constants.  

Figure 5

 

Comparison of cone density and spacing results calculated from automated MSCV-based cone identifications in five subjects with histology and published data (see Table 1). The data corresponding to five points (black arrows), one from each subject, are shown in Figure 6.

Figure 6

 

The retinal regions that were used to calculate cone spacing from five different subjects corresponding to the arrows in Figure 5B are shown. Images are ordered according to increasing cone spacing from (A) to (E) (i.e., listed in the order of arrows from bottom to top) with eccentricities of 5.52, 5.11, 5.58, 5.60, and 5.56 mm. Scale bar: 10 μm.

Both the histology and published data were within the range of the 95% prediction interval (PI) curves of the fit (Fig. 5B). Together these comparisons demonstrate that MSCV can be used to derive quick estimations of quantitative metrics such as cone density and spacing. 

Consistent Identification Results Between Confocal and Split Detection Images

Automated cone detections on confocal and split detection images were consistent with each other in five subjects, with average recall, precision, and F1 score of 87.1%, 88.9%, and 87.9% (Table 3; Fig. 7). These values are slightly smaller than the values reported in the earlier validation, but can be explained through close examination of the original images: additional bright circular structures were detected in the confocal image that did not have a corresponding cone-like object in the split detection images. In addition, the close proximity of cone photoreceptor cells reduces the overall contrast in the split detection images. 

Table 3

 

Evaluation of Identification Consistency Between Confocal and Split Detection Images

 

Figure 7

 

Comparison of cone identifications on confocal (A, C) and split detection (B, D) images from subject 9 at the eccentricity of 0.64 mm. Zooms of the white and black squares are shown (E–H). True identifications, green crosses; false positives, orange crosses; false negatives, yellow triangles, assuming that the automated confocal cone identifications are ground truth. Some false negatives (white circles) were caused by bright, circular structures that were only present in the confocal images but not present in the split detection images. These additional bright structures are possibly from the multiple reflections of a single cone in confocal images, or from rod photoreceptors. The yellow box in (A, B) shows where spacing and density was computed. Scale bars: 15 μm.

Automated identifications were also used to compute quantitative metrics based on a 55 × 55 μm ROI selected from the centers of each confocal-split detection image pair (yellow square in Fig. 7). The average relative errors of cone spacing and density were 1.2% and 10.0%, respectively (Table 4). 

Table 4

 

Comparison of Cone Spacing and Density Between Five Pairs of Confocal and Split Detection AOSLO Images in Table 3

 

Discussion

MSCV allows automated identification of cone photoreceptors in split detection AOSLO images. Its multiscale search makes it robust to variations in cone size with good results shown across eccentricities ranging from approximately 0.4 to 6.0 mm. MSCV-derived measurements of cone density and spacing are consistent with both previously published AO ophthalmoscopy and histology (Fig. 5). The overall robustness of the MSCV algorithm demonstrated on a full spectrum of split detection AO images showing variability in image quality suggests that it will be useful for establishing a larger normal database as well as for upcoming clinical applications in patients with retinal diseases. 

MSCV achieves high identification accuracy. As expected, image quality was the most important determining factor for identification accuracy. Poor image contrast caused MSCV to generate false positives and false negatives (Figs. 4D–F). We have suggested a new image quality assessment strategy based on anisotropy⁴⁸ that could be useful for evaluating split detection images with potential applications for filtering out low-quality AO images to ensure automated identification robustness. We envision that MSCV will significantly reduce the amount of time needed to assemble a more extensive AO-based database for photoreceptor metrics that includes both larger numbers of subjects across different age groups as well as data from a larger range of eccentricities, which will be particularly important for upcoming clinical trials that utilize AO as an outcome. Current AO-based normal databases are either lacking in sample size or in the range of eccentricities represented (Table 1), ostensibly due to the meticulous and tedious task of cone counting as well as the difficulty in identifying cones reliably in eccentric locations without the aid of split detection. On average, the time required to identify cones in each ROI was about 0.9 seconds, an order of magnitude faster than the current alternative based on manual identification (all speed measurements performed using a personal computer [Windows 7 PC with Intel quad-core i7-3770 3.4 GHz CPU; CPU release date April 2012] and 16 GB of random access memory). 

One other method for automated cone photoreceptor identification on split detection images has been proposed, based on the application of global smoothing to eliminate irregularities in illuminance, followed by local detection of image points containing intensity extremes.⁵² The main differences between this approach and ours is the scale-selection strategy in this paper to handle the issue that cone photoreceptors present with a substantial range of different sizes (Fig. 1), which requires variable levels of smoothing in different image regions, and the addition of a centering step which refines the positioning of the detection onto the actual cell center by automatically straddling bright and dark regions. Accurate centering of the identification is particularly important for spacing measurements. 

To the best of our knowledge, we have compared for the first time automated cone photoreceptor identification on simultaneously-acquired confocal and split detection images (Fig. 7). Systematic examination of differences between confocal and split detection images of the photoreceptor mosaic in an automated manner could be valuable for interpreting the status of cone photoreceptors in which cones appear “dark” or dim in the confocal modality.^{20,21,53–56} Automated identification of cones in the two different modalities reveal a number of discrepancies (white circles, Fig. 7) which are unrelated to the ability of the algorithm to detect cone-like structures but, rather, arise from differences in the appearance of cones in one modality versus the other. Some of these discrepancies can be explained by the presence of rods near the fovea⁵⁷ visible on the confocal images, which can mislead the cone identification algorithm. 

There are some important areas for future improvements to the MSCV algorithm, including the reduction of false positives and the recovery of false negatives. For the former, the integration of cone identification with cone segmentation algorithms developed for split detection images⁵⁸ could eliminate the problems of multiple identifications within a single-cone photoreceptor (orange circle in Fig. 4B). For the latter, machine learning could be employed to build cone photoreceptor classifiers to recover false negatives, as has been previously demonstrated for kidney cancer,⁵⁹ prostate cancer,⁶⁰ and lymph node identification.⁶¹ 

Conclusions

MSCV can quickly and automatically detect cone photoreceptors on split detection AOSLO images through multiscale image analysis. Validation results demonstrate that MSCV robustly identifies cone photoreceptors over a large range of eccentricities while still achieving high identification accuracy. There was good agreement in cone identifications between simultaneously acquired pairs of confocal and split detection images. Quantitative metrics derived from automated identifications were also comparable with histology and published data. These promising results can potentially accelerate the creation of larger AO-based normal databases for the rapid, noninvasive assessment of retinal pathology at the cellular level. 

Acknowledgments

The authors thank Catherine Cukras, Henry Wiley, Wadih Zein, Angel Garced, Meg Gordon, John Rowan, Patrick Lopez, Sharon Yin, Christina Appleman, Gloria Babilonia-Ayukawa, and Denise Cunningham for assistance with clinical procedures; Robert Cooper, Pavan Tiruveedhula, and Austin Roorda for sharing cone photoreceptor analysis software; Yusufu Sulai and Ethan Rossi for technical assistance with adaptive optics instrumentation; and Howard Metger for custom machining for adaptive optics instrumentation. 

Supported in part by the Intramural Research Program of the NIH, National Eye Institute, the Glaucoma Research Foundation Catalyst for a Cure Initiative, and NIH Grant U01 EY025477. 

Disclosure: J. Liu, None; H. Jung, None; A. Dubra, P; J. Tam, None 

Appendix

Multiscale Circular Voting

Let the intensity function of the split detection image after motion stabilization be defined as Display Formula\(I(x,y)\). Its multiscale image representation, Display Formula\(L\left( {x,y;t} \right)\), is generated through convolution with the Gaussian function Display Formula\(G\left( {x,y;t} \right).\)  

\begin{equation}\tag{A1}L\left( {x,y;t} \right) = G\left( {x,y;t} \right)*I(x,y)\end{equation}

The scale parameter Display Formula\(t\) is defined as Display Formula\(t \in \left[ {0.5, \cdots ,7.7} \right]\), where Display Formula\({t_{k + 1}} = {\rm{\gamma }}{t_k}\), with Display Formula\(k\) denoting the index of scale level. One important property of the multiscale image representation is that Display Formula\(L\left( {x,y;t} \right)\) has theoretically infinite derivatives.²⁸ In particular, derivatives and convolutions are commutable, which allows one to express derivatives of Display Formula\(L\left( {x,y;t} \right)\) as the convolution of derivatives of Display Formula\(G\left( {x,y;t} \right)\) with Display Formula\(I(x,y)\). Since Display Formula\(G\left( {x,y;t} \right)\) is infinitely differentiable, so is Display Formula\(L\left( {x,y;t} \right)\). As suggested by Mikolajczyk,⁴² the incremental constant Display Formula\({\rm{\gamma }}\) should be selected from Display Formula\(\left[ {1.1,1.4} \right]\). In this work, we use Display Formula\({\rm{\gamma }} = 1.2\), which was selected as a tradeoff between having a smooth transition between two consecutive scale levels (large values may cause aliasing artifacts) and between computation time (smaller values result in additional computations of Display Formula\(L\left( {x,y;t} \right)\)). The scale parameters are set based on prior knowledge of all possible radii (in pixels) of dark and bright cell regions in split detection images across the range of eccentricities that will be analyzed and for all possible fields of views.  

At each scale level Display Formula\(t\), we first compute the image gradient magnitude Display Formula\(\nabla L(x,y;t)\), and identify the points Display Formula\(\left( {{x_s},{y_s}} \right)\) for which Display Formula\(\nabla L({x_s},{y_s};t) \gt \alpha \).  

\begin{equation}\tag{A2}\nabla L({x_s},{y_s};t) = \left[ {{{\left( {{{\partial L} \over {\partial x}}} \right)}^2} + {{\left( {{{\partial L} \over {\partial y}}} \right)}^2}} \right]^{1/2} \approx \Biggl[ {{{\left( {{{L\left( {{x_{s + 1}},{y_s};t} \right) - L({x_{s - 1}},{y_s};t)} \over 2}} \right)}^2} + {{\left( {{{L\left( {{x_s},{y_{s + 1}};t} \right) - L({x_s},{y_{s - 1}};t)} \over 2}} \right)}^2}} \Biggr]^{1/2} \end{equation}

We set Display Formula\(\alpha = 10\) according to the original work of circular voting,³⁸ which we empirically determined to be effective for our datasets. Future work will reveal whether further optimization of this value leads to better results. A recursive approach was used to numerically compute Gaussian image smoothing (Equation A1) and its derivatives (Equation A2).⁶² For each selected point, Display Formula\(\left( {{x_s},{y_s}} \right)\), a voting area, A, is created around each of these (Fig. A1).  

\begin{equation}\tag{A3}A\left( {{x_s},{y_s},{r_{min}},{r_{max}},\Delta } \right) = \left\{ {\left( {{x_s} \pm r{\rm {cos}}\varphi ,{y_s} \pm r{\rm {sin}}\varphi } \right)|{r_{min}} \le r \le {r_{max}},\left. {\theta - \Delta \le \varphi \le \theta + \Delta } \right\}} \right.\end{equation}

 

Figure A1

 

Circular voting area \(A\), shown in gray, as defined by Equation A3.

Here, ± controls the detected region type: positive is used for bright cell regions, and negative for dark regions. Display Formula\(\theta \) is the voting direction initialized with Display Formula\(\mathop d\limits^ \to = - \nabla L({x_s},{y_s};t)/||\nabla L({x_s},{y_s};t)||\). Display Formula\(\Delta \) defines the angular search range of the voting area. The size of the voting area is constrained by Display Formula\({r_{min}} = 1\) and Display Formula\({r_{max}} = 2.5t\) pixels. Note that Display Formula\({r_{max}}\) utilizes the scale parameter t defined above, and as such Display Formula\({r_{max}}\), includes all possible radii of cell regions. This definition of Display Formula\({r_{max}}\) enables the MSCV algorithm to adapt to different AOSLO field of views without the need to specify separate parameter settings because it can automatically search for cone photoreceptors with varying sizes. This definition sets the upper bound of the voting area (Display Formula\({r_{max}}\)) such that the voting area is large enough to contain the cell boundary, but not so large that it contains too many points belonging to boundaries of neighboring cells (noting that bright and dark regions are approximately half the size of the cell and that distance is measured from the starting point toward the cell center as depicted in Figure 3A). The lower bound is the smallest possible region in pixel space (1 pixel). 

Figure A2

 

Radial symmetry leads to accumulation of MSCV responses in the center of circle-like regions. The small unfilled black circles represent image points with large gradient magnitudes defined in Equation A3 (e.g., cell region boundaries). Voting areas are cast toward the concave direction. Due to the radial symmetry, image points near the region center (black dots) are more likely to accumulate large MSCV response values, while the points near region boundary (gray dots) will yield less votes with smaller response values.

The MSCV response at an image point Display Formula\((x,y)\) contributed by a single selected point is computed as  

\begin{equation}\tag{A4}J\left( {x,y,{r_{min}},{r_{max}},\Delta ,{x_s},{y_s}} \right) = \mathop \sum \limits_{\left( {u,v} \right) \in A} \vert\vert\nabla L\left( {x + u,y + v;t} \right)\vert\vert{\hat {G}} \left( {u,v;t,\mathop d\limits^ \to ,A} \right),\end{equation}

Here, Display Formula\({\hat {G}} (u,v;t,\mathop d\limits^ \to ,A)\) is a Gaussian function with variance Display Formula\(t\), masked by the voting area Display Formula\(A\), and oriented along the voting direction Display Formula\(\mathop d\limits^ \to \). Note that Display Formula\((u,v)\) is an image point in the voting area Display Formula\(A\), which is used to compute the MSCV response at the point Display Formula\((x,y)\). Display Formula\({\hat {G}} \) is only valid on the voting area Display Formula\(A\), which is different from the Gaussian function Display Formula\(G\) in Equation A1 that is defined over the entire image domain. The MSCV response value is essentially the convolution between image gradient magnitudes and the Gaussian function Display Formula\({\hat {G}}\) over the voting area Display Formula\(A\). The radial geometry of cone photoreceptors enhances the MSCV response in the region centers as illustrated in Figure A2.  

MSCV is an iterative searching process. The MSCV response image Display Formula\(V\left( {x,y,{r_{min}},{r_{max}},{\Delta _i}} \right)\) is zeroed at each iteration, where Display Formula\(i \in \left[ {0,N} \right]\) is the iteration index and Display Formula\({\Delta _i}\) is the angular search range with Display Formula\({\Delta _N} \lt {\Delta _i} \lt {\Delta _0},{\Delta _0} = {\Delta _{max}},{\Delta _N} = 0\). Combining all voting response from different selected points, the overall MSCV response at each image point in i-th iteration is defined as  

\begin{equation}\tag{A5}V\left( {x,y,{r_{min}},{r_{max}},{\Delta _i}} \right) = \mathop \sum \limits_{({x_s},{y_s}) \in S} J(x,y,{r_{min}},{r_{max}},{\Delta _i},{x_s},{y_s})\end{equation}

After each iteration, the point with the maximum MSCV response is found in each voting area (white x's, Fig. 3).  

\begin{equation}\tag{A6}\left( {{x_m},{y_m}} \right) = \mathop {{\rm{argmax}}}\limits_{(x,y) \in A} V\left( {x,y,{r_{min}},{r_{max}},{\Delta _i}} \right)\end{equation}

The voting direction is next updated by the point Display Formula\(\left( {{x_m},{y_m}} \right)\),  

\begin{equation}\tag{A7}\theta \left( {x,y} \right) = {{({x_m} - x,{y_m} - y)} \over {\sqrt {{{({x_m} - x)}^2} + {{({y_m} - y)}^2}} }}\end{equation}

The MSCV response is computed using the current angular search range Display Formula\({\Delta _i}\), and is iteratively updated by reducing Display Formula\({\Delta _i}\) to zero.  

Combining region identifications at different scale levels produce a cluster of points for each cone photoreceptor. The point that contains the highest MSCV response is selected as the final identification of region center from the point cluster. 

Region Size Determination

The region size is determined through scale selection, using the scale-normalized Laplacian operator²⁸:  

\begin{equation}\tag{A8}\nabla _{norm}^2L\left( {x,y;t} \right) = t\left( {{{{\partial ^2}L} \over {\partial {x^2}}} + {{{\partial ^2}L} \over {\partial {y^2}}}} \right)\end{equation}

Its numerical computation can be simplified as the difference of multiscale image representations between two consecutive scales in terms of the heat equation, which can be derived as  

\begin{equation}\left( {{{{\partial ^2}L\left( {x,y;{t_k}} \right)} \over {\partial {x^2}}} + {{{\partial ^2}L\left( {x,y;{t_k}} \right)} \over {\partial {y^2}}}} \right) = {{\partial L\left( {x,y;{t_k}} \right)} \over {\partial t}}\end{equation}

 

\begin{equation}\tag{A9} \approx {{L\left( {x,y;{t_{k + 1}}} \right) - L\left( {x,y;{t_k}} \right)} \over {{t_{k + 1}} - {t_k}}} = {{L\left( {x,y;{\rm{\gamma }}{t_k}} \right) - L\left( {x,y;{t_k}} \right)} \over {{\rm{\gamma }}{t_k} - {t_k}}}\end{equation}

Therefore,  

\begin{equation}\tag{A10}\nabla _{norm}^2L\left( {x,y;{t_k}} \right) \approx {{L\left( {x,y;\gamma {t_k}} \right) - L\left( {x,y;{t_k}} \right)} \over {\gamma - 1}}\end{equation}

The optimal scale Display Formula\({\hat {t}}\) of each region is determined by searching for the most salient extrema:  

\begin{equation}\tag{A11}{\hat {t}}= {\rm{argmaxmi}}{{\rm{n}}_t}\nabla _{norm}^2L\left( {x,y;t} \right)\end{equation}

 

Centering

Let the coordinates and optimal scales of dual (bright and dark) region centers be Display Formula\(({x_d},{y_d},{{\hat {t}}_d})\) and Display Formula\(\left( {{x_b},{y_b},{{{\hat {t}}}_b}} \right)\). The location of the final cone identification, which is pushed closer to the region center with larger size, is calculated as  

\begin{equation}\tag{A12}x = {{{x_d}{{{\hat {t}}}_d} + {x_b}{{{\hat {t}}}_b}} \over {{{{\hat {t}}}_d} + {{{\hat {t}}}_b}}},y = {{{y_d}{{{\hat {t}}}_d} + {y_b}{{{\hat {t}}}_b}} \over {{{{\hat {t}}}_d} + {{{\hat {t}}}_b}}}\end{equation}