purpose. A methodology is proposed to build population-based average three dimensional (3-D) atlases or standards of the human cornea based on topographic data, along with variation maps. Also, methodologies for comparing populations or screening populations, based on these atlases are proposed.

methods. Topographies (Orbscan II; Bausch & Lomb, Rochester, NY) of 516 normal subjects were used. Methodology for the construction of a corneal atlas consisted of (1) data acquisition from both anterior and posterior corneal surfaces in the format of a 101 × 101 grid of *z* elevations evenly spaced (every 0.1 mm) along the *x* and *y* axes; (2) spatial normalization of the topographies on a unique average best-fit sphere to reduce the large variability in size and spatial location between corneas; (3) generation of the average 3-D model; and (4) statistics maps including average, median, and SD for each point of the grid.

results. To demonstrate the informative potential of this methodology, examples of atlases were generated. Numerical corneal atlases allow (1) characterization of a population, (2) comparison of two or more populations, (3) comparison of an individual with a reference population, and (4) screening of a population for the detection of specific corneal shape abnormalities, such as keratoconus or previous refractive surgery.

conclusions. The concept of a 3-D corneal atlas was developed. The proposed technique was meant to be simple, accurate, reliable, and robust and can be extended easily to any type of topographer capable of providing tridimensional corneal maps.

^{ 1 }

^{ 2 }generated average, standard deviation (SD) and count maps from the TMS-1 topographies of a group of 104 subjects. They also provided difference maps between average maps. Their goal was to compare superolateral versus superior phacoemulsification incisions. Buehren et al.

^{ 3 }facilitated the comparison between two populations by providing a Student’s

*t*-test map. Fam et al.

^{ 4 }also generated SD maps to evaluate the repeatability of the Orbscan pachymetry maps that they used to generate an average corneal pachymetry model. The repeatability was greatest in the center and decreased toward the periphery. The notion of alignment of the topographies before their averaging was raised by Buehren et al.

^{ 5 }who were concerned about eye movement between successive topographies. In their technique, a regression plane was used to correct for the tilt between two topographies, the best-fit sphere (BFS) apex was used to adjust for

*x*,

*y*,

*z*shifts, and a best-fit spherocylinder alignment was used to rectify the cyclodeviation error. In 2005, Grzybowski et al.

^{ 6 }addressed the usefulness of three different fitting zones (apex, global, and peripheral fit) for the alignment of pre- and post-LASIK topographies. They concluded that the choice of the fitting zone could influence the appearance of the central posterior elevation after LASIK. Dealing with right and left eyes also necessitates data conversion to account for the natural symmetry between both eyes with respect to the sagittal plane (enantiomorphism). The temporal cornea being steeper than the nasal cornea,

^{ 7 }simple averaging (without alignment of the corneas) results in an increased variability. Topuz et al.

^{ 8 }and Smolek and Klyce (Smolek MK et al.

*IOVS*2001;42:ARVO Abstract 2839) used the mirrored images of left eyes to convert them in right eyes before merging the data from both eyes.

^{ 1 }Topuz et al.,

^{ 8 }Fam et al.,

^{ 4 }and Smolek and Klyce groups (Smolek MK et al.

*IOVS*2001;42:ARVO Abstract 2839). However, alignment of the topographies before averaging was only used for topographies belonging to the same subject, and in none of the cases was the methodology for the development and optimization of these average models reported.

*x*,

*y*) of the 101 × 101 grid, the

*z*coordinate indicates the distance from the corneal surface to a reference plane perpendicular to the line of sight. Since the reference plane may vary from one topography to another, we made use of the anterior (or posterior) best-fit sphere (BFS) (i.e., the sphere that best adjusts to the anterior, or posterior, surface of the cornea in the sense of the least mean squares in a central adjustment zone of 10.0 mm in diameter). The BFS is characterized by its radius and center. Once the BFS was calculated from the elevation with respect to the reference plane, we used the elevation with respect to the BFS itself. We will explain in the next section how to compare two topographies using the BFS. The BFS also facilitates visualization of the corneal shape—namely, the highlighting of small elevation changes. Note that Orbscan provides its own BFS, but we preferred to compute the BFS ourselves, based on the raw elevation maps.

^{ 9 }

^{ 10 }This acoustic factor reduces all pachymetry values by 8% (Orbscan default acoustic factor of 0.92). In this study, we computed pachymetry according to the same method, using the acoustic correction factor.

_{ant/avg}of all radii obtained and the average coordinates (

*x*

_{ant/}

_{avg},

*y*

_{ant/}

_{avg},

*z*

_{ant/}

_{avg}) of all corresponding centers are then computed.

_{ant/avg}, as follows:

*x*,

*y*, and

*z*axes) of the surface points, to normalize the BFS radius to r

_{ant/avg}.

*x*

_{ant/}

_{avg},

*y*

_{ant/}

_{avg},

*z*

_{ant/}

_{avg}).

^{ 11 }This technique calculates a smooth surface passing by all realigned points and retrieves the surface elevation for the predefined grid points’ position.

*x*

_{ant/}

_{avg},

*y*

_{ant/}

_{avg},

*z*

_{ant/}

_{avg}) and the acoustic factor is taken into account, as just discussed.

*x*,

*y*) of the 101 × 101 grid, the average elevation and SD are calculated. As some of the points can be missing on the topographies, usually in the periphery (lashes, lid borders, and surface irregularities are usually responsible for this loss of information), statistics are computed only for the points for which at least 25% of the data are available.

^{ 12 }Large tear meniscus, surface irregularities, poor fixation by the subject, and/or incorrect alignment by the technician during acquisition are several possible causes of artifacts.

- Construction of an atlas for each population, as described herein.
- Computation of the elevation difference between the two atlases at each point (
*x*,*y*) of the 101 × 101 grid. - Computation of the appropriate
*P*statistics: For the comparison of two atlases, the Student’s*t*-test or the midrank Wilcoxon test^{ 13 }probability maps can be used, combined with the technique of Benjamini and Hochberg^{ 14 }to control the false-discovery rate. For the comparison of more than two atlases, this typical process can be generalized, and the computation of the probability can be replaced by an ANOVA or by an analysis of covariance, with the adjustment of the atlases values with one (or more) covariate(s) (age or time, for instance).

- Normalization of the tested topographic surface to fit the atlas average BFS.
- Computation of the point-based difference between the normalized tested topography and the atlas.
- Illustration with appropriate difference maps. Three difference illustration techniques are proposed: (1) The simple difference map is the computation of the point-based difference between the normalized map of the subject and the atlas average map. (2) The SD difference map is the computation of the point-based difference between the normalized map of the subject and the atlas average map, but where all differences that are less then three SD are set equal to zero. Thus, all points where the difference is greater than three SD from the average value will stand out. (3) The percentile map illustrates the subject corneal shape position with respect to the topographies composing the atlas. The percentile analysis requires, for each point of the 101 × 101 grid, the sorting of all
*z*elevation values of the topographies that were used in building the atlas. Formula 1 presents the percentile value perc as a function of the elevation*z*, whereas formula 2 gives the*z*elevation as a function of the percentile value:

*N*

_{ x,y }is the number of topographic maps covering the point (

*x*,

*y*). The parameter perc is the percentile between 1 and 100. A perc of 1 means that the elevation of the tested cornea, at that point, is at the level of the lowest

*z*elevation observed among the topographies used to build the atlas, whereas a perc of 100 means that, at that point, the elevation of the tested cornea is at the level of the highest

*z*elevation.

^{ 15 }based on shape indices derived from differences between the reference atlas and the tested subject, to automate the screening process.

_{ant/avg}of 7.94 mm and center coordinates (

*x*

_{ant/}

_{avg},

*y*

_{ant/}

_{avg},

*z*

_{ant/}

_{avg}) of (0.01, 0.01, −5.10) mm. The central corneal elevation was slightly positive (yellow-orange), with a mean elevation of 7.8 μm above the anterior BFS in the 2.0-mm diameter central zone). The SD was lower in the central region and increased toward the periphery (changing from cold to warm colors). For the posterior elevation map average BFS, we obtained r

_{post/avg}=6.59 mm and (

*x*

_{post/}

_{avg},

*y*

_{post/}

_{avg},

*z*

_{post/}

_{avg}) = (−0.02, −0.02, −4.37) mm. The posterior surface showed a more negative coefficient Q (increased prolateness), which resulted in an increased positive elevation centrally (23.7 μm above the posterior BFS in the 2.0-mm-diameter central zone). The SD map showed a slightly higher variability than for the anterior surface and an increased variability in the periphery (which could be attributed in part to a greater intersubject variability and a lower accuracy of the Orbscan for posterior surface and peripheral measurements.

^{ 16 }

^{ 17 }The average pachymetry map (Fig. 1A)showed a thinner central region. The thinnest point measured 547 μm and was slightly displaced temporally and inferiorly with respect to the topographic center (−0.2, −0.1). The variability was lower in the central region than in the periphery, a finding also reported by others.

^{ 4 }

*n*= 89 OD) and older (55–60 years,

*n*= 136 OD) subjects were studied. Overall, taken independently, both atlases looked quite similar, whether anterior elevation, posterior elevation, or pachymetry were considered. However, by subtracting the 20- to 25-year-old atlas from the 55- to 60-year-old atlas, the following observations were made. The anterior elevation difference map (Fig. 3 , top, third from left) revealed that the superior and inferior periphery was higher between 55 and 60 years than between 20 and 25 years, whereas the reverse was true of the nasal and temporal periphery. This corroborates the well-documented shift in astigmatism observed with age, from with the rule to against the rule.

^{ 18 }The posterior elevation difference map (Fig. 3 , center, third from left) revealed a slight temporal shift of the central bulge with age. Finally, the pachymetry difference map displayed a slight increase in the central pachymetry with age (Fig. 3 , bottom, 3rd from left).

*t*-test map, or mid-rank Wilcoxon test with the percentile information, with or without Benjamini’s correction, ANOVA).

^{ 19 }When subjects are compared, it has been shown that the major problem is excessive shape variability and that without alignment, the usefulness of the atlas is limited.

^{ 19 }

^{ 20 }

^{ 21 }Spatial normalization (e.g., translation, rotation, and scaling) was used and validated in brain imaging for shape analysis of brain tissues

^{ 19 }

^{ 20 }

^{ 22 }

^{ 23 }and for the construction of population-based lung shape average models.

^{ 24 }

^{ 25 }

^{ 26 }spatial normalization should be avoided, as the corneal refractive power would be altered proportionally to the stretching factor (scaling).

^{ 6 }or the periphery (e.g., 8.8- to 9.0-mm-diameter annulus). Again, this yields SD maps with low variations near the registration locus but much higher variations anywhere else.

Isotropic Scaling (%) | X,Y Translation (μm) | |||||
---|---|---|---|---|---|---|

Avg. ± SD | (Min; Median; Max) | Avg. ± SD | (Min; Median; Max) | |||

Anterior elevation | 0.00 ± 3.16 | (−11.06; −0.20; 12.63) | 18.37 ± 10.50 | (0.38; 16.27; 77.46) | ||

Posterior elevation | 0.00 ± 3.95 | (−11.25; 0.04; 14.83) | 48.84 ± 23.10 | (3.60; 46.93; 135.11) |

Isotropic Scaling (μm) Avg. ± SD (max) | X,Y Translation (μm) Avg. ± SD (max) | Scaling and Translation (μm) Avg. ± SD (max) | |
---|---|---|---|

Anterior elevation | 0.18 ± 0.01 (0.93) | 0.07 ± 0.01 (0.31) | 0.20 ± 0.01 (0.96) |

Posterior elevation | 0.83 ± 0.03 (2.39) | 0.20 ± 0.01 (0.60) | 0.85 ± 0.03 (2.50) |

Pachymetry | 3.87 ± 0.17 (24.67) | 0.11 ± 0.01 (0.51) | 2.26 ± 0.09 (12.53) |

**Figure 1.**

**Figure 1.**

**Figure 2.**

**Figure 2.**

**Figure 3.**

**Figure 3.**

**Figure 4.**

**Figure 4.**