## 1.

## Introduction

Corneal topography is a powerful clinical tool to evaluate corneal shape anomalies. Accurate methods and models are needed for the analysis of corneal topography and to evaluate their geometric properties,^{1} refractive errors, etc. Realistic corneal models are essential even for the reconstruction of the corneal surface from videokeratoscopy data.^{2, 3} Different parametric models of the corneal surface have been proposed to describe its elevation topography. Most of these models consider the sum of two terms: a regular basis surface (sphere, conicoid, biconic) plus an irregular component (polynomial expansion, splines, Fourier series, etc.).^{4, 5, 6} These parametric models can adjust and describe the topography of normal corneas with reasonable or even high accuracy. However, in cases of severe irregularities or deformations (keratoconus, postsurgical corneas, etc.), these models may fail to describe the real topography.^{7} Nevertheless, difficulties to fit models to abnormal corneas can be realized as a potential tool to classify between normal or pathological corneas. Some parameters of corneal models can be used as descriptors for automatic classification procedures. In this sense, several successful and clinically relevant procedures have been developed to detect incipient keratoconus or classify corneas.^{8, 9, 10, 11, 12}

In particular, postsurgical corneas may present discontinuities (either in the surface or in its first or second derivatives), and hence realistic models should be able to show these discontinuities as well. In refractive surgery, ablated corneas present different zones with different geometrical and refractive properties. Standard treatments apply spherical-cylindrical ablation profiles, whereas modern customized treatments use a wider variety of ablation profiles, including asymmetrical ones, based on optical and anatomical characteristics of the patient’s eye. For instance, myopic LASIK correction mainly works in the central optical zone (OZ), but laser ablation is also applied to the transition zone (TZ), and only a peripheral corneal ring (P) remains untouched. The goal of the central ablation pattern is to reduce the curvature in the OZ, while the ablation at TZ smooths the curvature changes (transition), thus avoiding discontinuities. Since the ablation patterns at OZ, TZ, and P are different, significant differences should be expected in the resulting postsurgical topography among these three zones. Typical sizes of these zones are as follows: OZ—about $6\text{-}\mathrm{mm}$ diameter; TZ— $6\text{-}\mathrm{mm}$ inner and $9\text{-}\mathrm{mm}$ outer diameter; and P—from about $9\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ to the limbus $(11\u201312\phantom{\rule{0.3em}{0ex}}\mathrm{mm})$ .

Our departure hypothesis is that ablated corneas must somehow exhibit these well-differentiated zones; hence, the analysis of corneal topography should be made independently for each zone. Therefore, the aim of this work is to develop an automatic method to build custom multizone models of postsurgical corneas and then to apply this procedure to analyze and compare the outcome of three different, standard, and custom LASIK treatments. To build the multizone model, we analyze the corneal elevation topography in two stages:

1. Surface segmentation into a predefined number of zones (three in myopic LASIK) is performed, based on a set of point descriptors and a clustering algorithm. In this way, each point of the surface is assigned to a particular zone based on the descriptor statistics.

2. Once segmented, a standard model (nonrevolution conicoid plus Zernike polynomial) is applied to each zone (optical, transition, and periphery) independently.

^{13}It must be noticed that computing the difference between pre- and postsurgical elevation maps directly does not generally allow reliable ablation profiles to be obtained,

^{14}since they are referred to planes or surfaces that generally do not match.

The proposed multizone model was applied to analyze and compare the outcome of three different systems of LASIK surgery. One of them is for standard treatment and the other two systems are for custom myopic LASIK. For this purpose, corneal models were obtained from both preoperative and post-LASIK elevation topographies for each patient in three groups corresponding to the different custom and standard treatments. Subsequent morphological and optical computations permitted us to obtain relevant parameters of the optical zone such as the shape and center of the ablation, curvatures, conic constants, optical quality, etc., and then to compare programmed versus implemented values. We also studied analogies and differences between the outcomes of the different treatments.

## 2.

## Material and Methods

The procedure to build the multizone model of ablated corneas is depicted in Fig. 1 . The main stages are

1. Fitting a standard monozone model to the input elevation topography data.

2. Computing a set of local descriptors able to discriminate among the different zones. For each point of the corneal surface, these descriptors are the curvature, the fit error, and the distance to the center.

3. Performing a segmentation of the corneal surface into three zones (OZ, TZ, and P) using a clustering algorithm based on the above descriptor set.

4. Fitting again the standard model, but now to each of the three zones independently.

5. The final model will be the union of the three zones.

## 2.1.

### Standard Model

In most models, the corneal surface $z=S(x,y)$ can be expressed as the sum of two terms:

where $b(x,y)$ is a regular basis surface describing the overall shape of the cornea, and $r(x,y)$ accounts for irregularities and departures from the basis surface. If the basis surface $b(x,y)$ is an accurate model of the average topography of the cornea, then $r(x,y)$ will tend to be low, at least for normal corneas, and hence the analysis of the corneal surface $b(x,y)$ will provide meaningful overall properties. The residual, $r(x,y)=S(x,y)-b(x,y)$ , is usually fit to some polynomial expansion. The most widely used one is the Zernike polynomial^{5}given bywhere ${P}_{k}$ is the $k\text{th}$ Zernike polynomial, and ${c}_{k}$ is its corresponding coefficient. Here we follow the OSA standards for reporting the optical aberrations of eyes, and consider them to the seventh order,

^{15, 16}that is, 36 terms (including 0th order, piston).

We consider a general conicoid as our basis surface
$b(x,y)$
, as it provided highly satisfactory results to describe the overall shape of normal corneas:^{2}

## Eq. 3

$${a}_{11}{x}^{2}+{a}_{22}{y}^{2}+{a}_{33}{z}^{2}+{a}_{12}xy+{a}_{13}xz+{a}_{23}yz+{a}_{1}x+{a}_{2}y+{a}_{3}z+{a}_{0}=0.$$In what follows, the input data are Orbscan II elevation topography maps. A more detailed description of the model and least-squares fit can be found in Ref. 2. This model as well as all subsequent computations were implemented using Matlab (The Mathworks Inc., Natick, MA).

## 2.2.

### Local Descriptors

Due to laser ablation, post-LASIK corneas can have different zones, characterized by different properties, and eventually they can be separated by some discontinuity. In these cases, the standard model could fail to fit the topography around discontinuities.^{7} The idea here is to overcome this drawback and use all the available information to identify, point by point, the different zones according to homogeneous features.

Corneal descriptors will be those parameters, defined for every point of the surface, which can be used to classify or assign each point to a given zone. Ideally, descriptors must be homogeneous within each zone but must differ between zones. In addition, the number of descriptors must be the minimum able to guarantee an efficient discrimination. In our case, we found that this minimum number was three. To be efficient, these parameters must be independent from each other, and hence describe totally different properties. For each point, the chosen descriptors were distance to the center of the cornea $\left(D\right)$ , curvature $\left(C\right)$ , and fit error $\left(E\right)$ .

## 2.2.1.

#### Distance to the center, $D$

Here we use the a priori knowledge of the characteristic geometry of myopic LASIK treatments, with a central, OZ, and two successive concentric rings, TZ and P. Thus, it is evident that for low values of $D$ , there is a high probability that the point belongs to the OZ; for intermediate values, the highest probability is to belong to the TZ, and to P for the highest values. Therefore, this descriptor can help to discriminate among points of the different zones.

## 2.2.2.

#### Local curvature, $C$

The objective of LASIK treatments is correcting the refractive error by modifying the corneal curvature. It is known that the corneal curvature is not totally homogeneous and changes from point to point even in normal corneas. Nevertheless, it seems reasonable to assume that after LASIK the resulting curvature will probably be more homogeneous within a zone than between different zones. In particular, the optical zone is expected to have a more homogeneous curvature, whereas transitional or peripheral zones may display a higher variability. Our descriptor will be the Gauss curvature,
$C$
, which is the product of the principal curvatures,
${k}_{1}$
and
${k}_{2}$
, defined as the maximum and minimum curvatures among all orientations. Mathematically, the Gauss curvature was computed for each corneal point from coefficients of the first and second fundamental forms^{17} of the corneal surface according to the following expression:

^{18}$K=1000(1-1.3375)\sqrt{C}$ , assuming 1.3375 for the refractive index of the cornea.

## 2.2.3.

#### Fit error, $E$

The third descriptor measures the departure from the adjusted standard model to the actual topography. As a way to enhance the role of this descriptor, here we directly consider the modulus of the residual
$E=\mid r(x,y)\mid $
obtained after fitting the basis surface (ellipsoid) alone. Roughly speaking,
$E(x,y)$
is a measure of the local irregularity. Similarly to the other descriptors, one may expect that the magnitude of the fit error should be more homogeneous within each zone but change substantially between zones. The magnitude of this corneal descriptor will depend on the basis surface chosen (sphere, ellipsoid, etc.), but as an *irregularity* descriptor its discrimination capability may be relatively invariant against changes to the basis surface.

This choice of descriptors is probably not unique, but here we have tried to apply the following criteria: The descriptors must be meaningful physical magnitudes, independent, easy to compute, and able to discriminate between zones. Note that the last requirement is especially important and requires empirical validation, which is one of the primary goals of this study.

These descriptors are used as the three components of a feature vector ( $D$ , $C$ , $E$ ) at each point of the corneal topography. As they represent different magnitudes and can have totally different ranges, it is necessary to normalize each descriptor to its mean value in order to have dimensionless variables with the same weight in the segmentation. After normalization, we can also give different weights to enhance or reduce the relative influence of each descriptor. We tried different weights, but the results did not change substantially. We finally applied the weights (1, 2, 1) to take into account that the Gauss curvature is probably the most important feature.

## 2.3.

### Segmentation

Once the
$n=3$
descriptors are calculated for every point of the corneal surface, we build an
$n$
-column matrix, where each row corresponds to one corneal point, and each column to one descriptor value. Then segmentation is carried out by applying a
$k$
-means clustering algorithm. This type of algorithm tries to find
$k$
clusters in the three-dimensional histogram of the
$n$
-column matrix, where
$k$
is the number of classes. For convenience, in our implementation we consider
$k=4$
: namely, the three zones OZ, TZ, and P, plus an additional class corresponding to the outer area formed by those points having null value (no topography data available). Among the different possible metrics, we have selected the Euclidean distance between points in the 3D histogram. The implementation used the Matlab
$k$
-*means* function. The result of the segmentation is that each point in the topography is assigned to one of the three zones. In other words, we obtain three binary maps or masks,
${M}_{OZ}$
,
${M}_{OT}$
, and
${M}_{P}$
, corresponding to the optical zone, transition region, and periphery, respectively. Each binary map takes value 1 within its region and 0 elsewhere. In this way, the ablated cornea is divided into three well-defined spatial regions, which we can analyze independently. The
$k$
-means algorithm does not guarantee that the three zones are topologically correct and connected a priori, and the output could contain some misclassified isolated points. However, we never obtained disconnected areas or isolated points with this set of descriptors, so that further segmentation refining was unnecessary.

## 2.4.

### Zonal Analysis

The final stage is to fit the standard model of Eq. 1 to each zone, defined by the corresponding mask. As a result, each corneal zone is characterized by three functions; the mask $M$ providing the area, the basis ellipsoid surface $b(x,y)$ , and the residual $r(x,y)$ . These three functions permit one to perform a more detailed analysis of each zone. Here, we will focus mainly on the optical zone. The basis function provides curvature radii $R$ and conic constants $Q$ along the $x$ - and $y$ -axes, as well as the position and orientation of the optical axis; the residual analysis yields the Zernike coefficients ${c}_{k}$ . Apart from the standard analysis, the present multizone approach adds new useful information by simple standard morphological analysis of the masks. In the particular case of the optical zone, we can simply perform a least-squares fit of the perimeter of the optical zone to a free oriented ellipse in the 2D space with semiaxes ${a}_{OZ}$ and ${b}_{OZ}$ . This fit to the mask ${M}_{OZ}$ provides three interesting parameters:

1. The

*optical zone diameter*can be computed either as the diameter of the best-fit circumference to the perimeter of the optical zone or as $\sqrt{{a}_{OZ}^{2}+{b}_{OZ}^{2}}$ .2. The

*eccentricity*is defined by ${\epsilon}^{2}=1-{b}_{OZ}^{2}\u2215{a}_{OZ}^{2}$ . When $\epsilon =0$ , we have a circumference; as it departs from that value, we have a progressively elongated ellipse.3. The

*OZ center*is directly obtained as the center ${x}_{0}$ , ${y}_{0}$ of the adjusted ellipse (alternatively, it can be computed as the centroid of the optical zone). The center of coordinates represents displacements from the center of the topography map. In our case, this last term corresponds to the vertex normal, i.e., the intersection of the keratometric axis with the anterior corneal surface. Thus, ${x}_{0}$ , ${y}_{0}$ represent decentrations of the ablation pattern with respect to the vertex normal.

## 2.5.

### Comparative Study

The above methods were applied to a comparative study of three different (standard and customized) myopic LASIK techniques. A total of 31 eyes (25 patients) were analyzed. All subjects (15 female, 10 male), aged between 21 and $64\phantom{\rule{0.3em}{0ex}}\text{years}$ (mean $35\pm 10\phantom{\rule{0.3em}{0ex}}\text{years}$ ), were regular patients in the RealVision clinic (Madrid) and had their treatments by the same surgeon. All patients were myopic with spherical refractions in the range $-0.25\phantom{\rule{0.3em}{0ex}}\mathrm{D}\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}-7.5\phantom{\rule{0.3em}{0ex}}\mathrm{D}$ (mean $-4.3\pm 2.2\phantom{\rule{0.3em}{0ex}}\mathrm{D}$ ). Cylinder refraction ranged from $0\phantom{\rule{0.3em}{0ex}}\mathrm{D}\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}-3\phantom{\rule{0.3em}{0ex}}\mathrm{D}$ (mean $-0.9\pm 0.6$ ). The three groups were:

• Group 1: Ten eyes had standard surgery with a PlanoScan™ 2000 (Technolas®217A excimer laser, Bausch&Lomb, Salt Lake City, UT) system. The average initial refraction was $-3.5\pm 2.4\phantom{\rule{0.3em}{0ex}}\mathrm{D}$ (sphere) and $-1.2\pm 0.8\phantom{\rule{0.3em}{0ex}}\mathrm{D}$ (cylinder).

• Group 2: Thirteen eyes had custom Zyoptix (Bausch&Lomb) treatment. Average refraction: $-4.6\pm 2.3\phantom{\rule{0.3em}{0ex}}\mathrm{D}$ (sphere) and $-0.9\pm 0.5\phantom{\rule{0.3em}{0ex}}\mathrm{D}$ (cylinder).

• Group 3: Eight eyes had custom treatment, with an Allegretto (WaveLight Laser Technologie, Erlangen, Germany) system. Average refraction: $-4.8\pm 2.1\phantom{\rule{0.3em}{0ex}}\mathrm{D}$ (sphere) and $-0.6\pm 0.6\phantom{\rule{0.3em}{0ex}}\mathrm{D}$ (cylinder).

For each eye, three preoperative and three postoperative (three months after surgery) corneal elevation topographies were taken using an ORBSCAN II apparatus. In each case, the average pre-LASIK and post-LASIK topographies were computed.^{2, 19, 20} A multizone model was obtained for each postsurgical cornea, and the analysis described above was applied.

## 3.

## Results

The results are organized in two parts. First, we present the complete analysis for one eye to illustrate the type of information provided by our method, whereas the second part summarizes the results for the three groups of patients.

## 3.1.

### One Example

The case analyzed here corresponds to the right eye of patient AFD, female, $34\phantom{\rule{0.3em}{0ex}}\text{years}$ old. She was treated to correct a refractive error of $-4\phantom{\rule{0.3em}{0ex}}\mathrm{D}$ sphere and $-0.50\phantom{\rule{0.3em}{0ex}}\mathrm{D}$ cylinder at 150°, with custom Allegretto LASIK (group 3) system. Figure 2 compares pre- (left) and postoperative (right) corneal elevation topographies, with respect to their best-fit spheres. We can see that the topography has totally changed with the surgery, and hence it is difficult to make a direct comparison between them.

To develop the multizone model, we first obtain the three descriptors, shown in Fig. 3 before normalization, with a negative gray scale (white corresponds to lowest values and black to highest values). The left panel represents $D$ , the distance to the center; the scale is in millimeters. This descriptor is constant for all topographies, as this distance is independent from them. As expected, it presents a quite regular and smooth pattern. The central panel corresponds to the Gauss curvature, $C$ , in diopters, and the right panel to the resulting fit error, $E$ , in micrometers. There is a high correlation between, $C$ and $E$ despite their different physical meaning and magnitudes (diopters and microns, respectively). In both cases, the spatial distribution of the different zones appears clearly represented. Note that this clear distinction can not be made directly on the elevation topography (Fig. 2, right panel). Interestingly, the Gauss curvature shows a somewhat granular appearance, resembling laser spot imprints. Such a granular appearance could be observed in most cases.

Each topography map consists of about 10,000 data points, and hence we have the same number of points per descriptor. Figure 4 shows the data points in the 3D descriptor space, which is the main input to the clustering algorithm. Here the descriptors were normalized to their mean values, and $C$ (Gauss curvature) was weighted by a factor of two. The average descriptor values are listed in Table 1 (upper rows) for the three zones, whereas the solid symbols in Fig. 4 represent normalized average values on each zone. Points of the transition zone are characterized by intermediate distances $D$ , somewhat higher errors $E$ , and a considerably higher curvature $C$ . Thus, TZ points tend to lie on the back, in the upper part, and in the middle between OZ (left) and P (right) points. The most distinctive feature of peripheral points is their high distances, so that they lie on the right side; descriptors $E$ and $C$ show a wider variability so that P points spread more along these axes. Conversely, distances are much lower in the optical zone, where the distribution of points is more regular. Average values in Table 1 clearly show these trends. The “monozone” column displays the global average of the three descriptors, which is close but not equal to the average over zones.

## Table 1

Main results of the analysis of corneal topographies for subject AFD. Length units are in millimeters, RMS (both Zernike and errors) are in micrometers, and angles in degrees.

Pre-LASIK | Monozone | OZ | TZ | P | |
---|---|---|---|---|---|

Average descriptors | |||||

D (mm) | 3.50 | 1.86 | 3.59 | 4.76 | |

C (diopters) | 42.05 | 38.46 | 46.42 | 38.97 | |

E $\left(\mu \mathrm{m}\right)$ | 0.11 | 0.10 | 0.12 | 0.10 | |

Zone morphology | |||||

Diameter | 5.46 | 8.66 | |||

Eccentricity | 0.25 | 0.30 | |||

${x}_{0}$ , ${y}_{0}$ (center) | $-0.17$ , $-0.44$ | $-0.01$ , 0.27 | |||

Surface topography | |||||

${R}_{x}$ | 7.63 | 8.61 | 9.37 | 8.41 | 6.51 |

${R}_{y}$ | 7.34 | 8.37 | 9.08 | 8.11 | 6.32 |

${Q}_{x}$ | $-0.36$ | 0.41 | 1.05 | 0.12 | $-0.91$ |

${Q}_{y}$ | $-0.39$ | 0.37 | 0.99 | 0.08 | $-0.91$ |

${x}_{a}$ , ${y}_{a}$ (apex) | $-0.04$ , $-0.59$ | 0.46, 0.33 | 1.13, 0.15 | ||

$\alpha $ $\beta $ (axis)1 | 0.3°, 4.5° | $-2.1\xb0$ , $-2.6\xb0$ | $-7.1\xb0$ , $-1.0\xb0$ | ||

RMS Zernikes | 11.50 | 13.11 | 1.47 | ||

RMS fit error | |||||

2.43 | 2.91 | 0.33 | 1.14 | 1.79 |

Figure 5 shows the result of the clustering algorithm. Each point in the topography map is assigned to one of the three possible zones, and the three resulting masks correspond to the optical zone, transition zone, and periphery. The outer white area corresponds to the fourth class, namely points with no elevation data.

After segmentation, we can complete the procedure depicted in Fig. 1 to obtain different subsets of results. First, we obtained the morphological data of the different zones in Fig. 5. These are listed in the three “zone morphology” rows of Table 1. We can see that the computed optical zone is not circular, but has some elongation (0.25 eccentricity), with an average diameter of
$5.46\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
. The TZ has a similar elongation (0.3), and an external diameter of
$8.66\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
. The center of the OZ (and possibly the ablation center) is displaced from the vertex normal
$0.44$
(nasal) and
$0.17\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
upwards. Here we want to remark that these resulting zones, OZ, TZ, P, are “effective” areas,^{21, 22, 23} defined according to the statistical distribution of the three descriptors. This means that the optical zone is an area where the curvature, the fit error, and the distance to the center have relatively homogeneous values, and the same applies to TZ and P. Conversely, points in OZ must have significantly different descriptor values than points in TZ. Therefore, as the boundary between two zones is established in statistical terms, the resulting zones can be considered as functionally homogeneous or effective, and hence they may differ from the planned nominal areas.

The next step was to fit our standard topography model to the different zones. The resulting parameters describing the surface topography are also listed in Table 1. The values obtained from the (monozone) analysis of the presurgical topography are given in the first column, and the same standard analysis for the postsurgical cornea is given in the second column. Within each column, differences between
${R}_{x}$
and
${R}_{y}$
correspond to corneal astigmatism. It is clear that the surgery increases the radii, mainly in the optical zone (
${R}_{x}$
from
$7.63\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
up to
$9.37\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
), thus reducing the power of the cornea, as expected. The conic constants show the typical well-known trend in LASIK surgery.^{13, 24} Normal corneas have negative Q values: in this case
$-0.36$
and
$-0.39$
, which are close to average values.^{2} However, the sign becomes positive after LASIK surgery, which is patent even for the monozone analysis
$(+0.41,+0.37)$
. Several studies have demonstrated that this causes the increase in spherical aberration after LASIK surgery.^{25, 26, 27}. Our multizone analysis reveals that the actual increase in the Q values in the optical zone is even greater (about
$+1$
). This will produce a strong increase in spherical aberration. On the other hand, the peripheral zone seems to keep negative Q values, as it was untouched by the laser. The apex and orientation of the optical axis also show significant changes after surgery, and again, the changes are more pronounced when we consider the optical zone alone. Before surgery, the orientation of the optical axis is natural (0.3° downwards, 4.5° nasal), but after LASIK these angles are reversed. The resulting OZ optical axis is tilted upwards. 7.1°. A further analysis of the shape of the surface is given by the Zernike coefficients, which account for irregularities and departures from the basis ellipsoid geometry. The monozone analysis shows around a 15% increase in the total Zernike RMS (from
$11.50\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}13.11\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
), and hence in the irregularity. The optical zone, however, seems to follow the ellipsoid geometry much better (RMS of
$1.47\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
). However, we cannot directly compare the RMS values of OZ with those of the monozone model, which considers the whole corneal topography.

Finally, the reconstruction of the corneal topography, as the union of the models obtained for the three zones (Fig. 6 , left panel), shows a high fidelity with the original data (right panel). The final RMS fit error obtained for all cases and zones is listed in the bottom row of Table 1. The RMS fit error is low in all cases (between 2 and $3\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ in pre-LASIK and post-LASIK monozone models). Nevertheless, the multizone model provides lower errors, especially for OZ $\left(0.33\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}\right)$ .

## 3.2.

### Comparative Study

The goodness of fit for the complete set of patients and models is summarized in Fig. 7 in terms of the RMS fit error. Error bars (standard deviation) represent subject variability. From this figure, it is patent that the Zyoptix group had lower fit error before surgery; on the contrary, the Allegretto group showed higher errors; and the PlanoScan group represents an intermediate case. Roughly speaking, a lower fit error would mean a more regular topography, and vice versa. After surgery, we find basically the same relationship; i.e., no significant differences in postsurgical changes among groups. There is a significant increase in the fitting error after LASIK for the monozone model (global average passes from $1.43\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}2.47\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ ). The optical zone, however, shows a much lower error (average $0.28\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ ), whereas TZ, P, and the whole multizone model (average $1.24\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ ) keep the error similar to that of the pre-LASIK level. These results are essentially the same obtained above for subject AFD.

Figure 8
compares the estimated versus programmed diameter of the optical zone. In all cases, the estimated, or effective,^{21} diameter is lower than the nominal programmed value. The average values were
$6.57\pm 0.27\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
and
$5.50\pm 0.40\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
, respectively, so that the difference is slightly above
$1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
(
$1.05\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
,
$1.15\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
, and
$0.97\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
for PlanoScan, Zyoptix, and Allegretto, respectively). As we said before, this difference is explained by the different definition of programmed and effective OZ. We cannot determine whether this effect is due to surgery or is a bias of our segmentation algorithm. Both factors probably contribute to reduce the effective diameter.

The eccentricity, or elongation, of the effective optical zone is $0.29\pm 0.13$ for PlanoScan, $0.24\pm 0.13$ for Zyoptix, and $0.31\pm 0.17$ for Allegretto (global average $0.27\pm 0.14$ ). The average elongation was nearly horizontal, as in the example of Fig. 5, but the orientation showed a large intersubject variability: $3\xb0\pm 44\xb0$ . Among groups, Zyoptix displays a slightly lower elongation, i.e., more circular OZ. Conversely, the Allegretto treatment seems to produce slightly larger elongations. The coordinates of the ablation center with respect to the vertex normal (origin of the corneal topography map) are given in Fig. 9 . Average decentrations are low for the three groups: ${x}_{0}=0.00\pm 0.06\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , ${y}_{0}=0.00\pm 0.09$ for PlanoScan; $0.01\pm 0.09\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , $0.02\pm 0.10$ for Zyoptix; $0.07\pm 0.24\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , $0.10\pm 0.10\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , for Allegretto; global average, ${x}_{0}=0.02\pm 0.13\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , ${y}_{0}=0.04\pm 0.10$ . Therefore, ablation profiles seem well centered and with quite reasonable standard deviation values, especially for the PlanoScan group. Only four cases (three Allegretto, one Zyoptix) showed decentrations greater than $0.25\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ . A possible explanation for these small differences found among groups might be related to the degree of the surgeon’s training with the different systems.

The analysis of the corneal topography model provides the curvature radii, conic constants, optical axis, and Zernike coefficients. In myopic LASIK, the curvature radii increase after surgery (see Table 2 ). As in the example, post-LASIK average values also show higher values in the optical zone, as compared to the monozone model; the transition zone shows a smaller increase with intermediate values, whereas the periphery shows lower radii. If we compare groups, both PlanoScan and Zyoptix present similar outcomes, but Allegretto shows higher increments. Interestingly, for the last group, the curvature radii do increase even in the periphery.

## Table 2

Average curvature radii (mm) along horizontal (H) and vertical (V) directions for the three groups of patients. Pre-LASIK and Post-LASIK correspond to monozone models and OZ, TZ, and P correspond to the three zones after LASIK.

Pre-LASIK | Post-LASIK | OZ | TZ | P | ||
---|---|---|---|---|---|---|

PlanoScan | HV | $7.66\pm 0.27$ $7.39\pm 0.21$ | $8.31\pm 0.24$ $8.06\pm 0.29$ | $8.46\pm 0.38$ $8.29\pm 0.41$ | $8.11\pm 0.24$ $7.91\pm 0.30$ | $6.98\pm 3.32$ $6.77\pm 3.20$ |

Zyoptix | HV | $7.65\pm 0.10$ $7.45\pm 0.16$ | $8.16\pm 0.33$ $7.96\pm 0.30$ | $8.48\pm 0.50$ $8.27\pm 0.39$ | $8.05\pm 0.31$ $7.82\pm 0.28$ | $7.18\pm 1.57$ $6.98\pm 1.52$ |

Allegretto | HV | $7.85\pm 0.30$ $7.59\pm 0.28$ | $8.98\pm 0.47$ $8.75\pm 0.42$ | $9.26\pm 0.70$ $9.12\pm 0.59$ | $8.68\pm 0.38$ $8.52\pm 0.44$ | $8.23\pm 0.84$ $7.94\pm 0.77$ |

Average | HV | $7.70\pm 0.24$ $7.47\pm 0.22$ | $8.42\pm 0.48$ $8.20\pm 0.46$ | $8.68\pm 0.62$ $8.50\pm 0.57$ | $8.23\pm 0.41$ $8.03\pm 0.44$ | $7.39\pm 2.17$ $7.16\pm 2.09$ |

The resulting average conic constants are listed in Table 3
. In this case, only values for
${Q}_{x}$
are given, because
${Q}_{x}$
and
${Q}_{y}$
are highly correlated in the ellipsoid model.^{2} Average values show the same trends as in the above example. The initially negative conic constants
$(-0.41)$
become positive after surgery. This effect is more marked in the optical zone (average
$+0.29$
). Among groups, Zyoptix provides lower Q values, but positive (average
$+0.14$
in OZ), whereas Allegretto produces higher values (average
$+0.55$
in OZ). As we will discuss later, passing from negative to positive Q values after surgery will notably increase spherical aberration.

## Table 3

Conic constant Qx (horizontal) for the three groups of patients.

Pre-LASIK | Post-LASIK | OZ | TZ | P | |
---|---|---|---|---|---|

PlanoScan | $-0.41\pm 0.17$ | $0.04\pm 0.40$ | $0.27\pm 0.62$ | $-0.14\pm 0.38$ | $-0.45\pm 0.92$ |

Zyoptix | $-0.41\pm 0.10$ | $-0.08\pm 0.33$ | $0.14\pm 0.61$ | $-0.18\pm 0.28$ | $-0.68\pm 0.83$ |

Allegretto | $-0.40\pm 0.07$ | $0.49\pm 0.22$ | $0.55\pm 0.31$ | $0.19\pm 0.16$ | $-0.08\pm 0.40$ |

Average | $-0.41\pm 0.12$ | $0.11\pm 0.40$ | $0.29\pm 0.56$ | $-0.07\pm 0.33$ | $-0.45\pm 0.79$ |

The average orientation (angles $\alpha $ and $\beta $ with respect to the $x$ and $y$ -axes, respectively) and coordinates of the corneal intercept, apex, of the corneal optical axis are listed in Table 4 . In ellipsoid model 2, the optical axis is defined by the orientation and position of the ellipsoid axis along $z$ . It is patent that the surgery modifies both the position and orientation. The angles $\alpha $ and $\beta $ show much higher standard deviations than averages, thus suggesting a wide intersubject variability. Such variability increases strongly after surgery (typical values above 10°), indicating low predictability. The apex coordinates also change after surgery, typically by several tenths of a millimeter. Again, standard deviations are higher than averages, so that there seems to be a strong random component.

## Table 4

Average orientation of the optical axis and apex position for the three groups of patients. For the pre-LASIK case, the data correspond to the global monozone model, whereas we only consider the optical zone for post-LASIK corneas. Angles are positive for downside and nasal orientations. Coordinates are in millimeters with respect to the vertex normal.

Pre-LASIK | OZ Post-LASIK | ||
---|---|---|---|

PlanoScan | $\alpha $ $\beta $ | $-1.09\xb0\pm 5.49\xb0$ $-0.01\xb0\pm 3.56\xb0$ | $1.71\xb0\pm 3.56\xb0$ $-3.94\xb0\pm 10.01\xb0$ |

${x}_{0}$ ${y}_{0}$ | $-0.09\pm 0.75$ $-0.01\pm 0.46$ | $0.14\pm 0.53$ $-0.54\pm 1.38$ | |

Zyoptix | $\alpha $ $\beta $ | $1.04\xb0\pm 2.19\xb0$ $0.20\xb0\pm 2.24\xb0$ | $6.87\xb0\pm 13.31\xb0$ $-0.64\xb0\pm 9.17\xb0$ |

${x}_{0}$ ${y}_{0}$ | $-0.09\pm 0.32$ $0.02\pm 0.30$ | $-0.77\pm 1.27$ $0.03\pm 1.29$ | |

Allegretto | $\alpha $ $\beta $ | $1.17\xb0\pm 1.39\xb0$ $-1.92\xb0\pm 2.47\xb0$ | $-2.44\xb0\pm 16.76\xb0$ $-0.58\xb0\pm 12.69\xb0$ |

${x}_{0}$ ${y}_{0}$ | $0.01\pm 0.26$ $-0.25\pm 0.33$ | $0.63\pm 1.70$ $0.30\pm 1.00$ | |

Average | $\alpha $ $\beta $ | $0.13\xb0\pm 3.99\xb0$ $-0.39\xb0\pm 2.99\xb0$ | $2.80\xb0\pm 12.45\xb0$ $-1.69\xb0\pm 10.19\xb0$ |

${x}_{0}$ ${y}_{0}$ | $-0.06\pm 0.53$ $-0.06\pm 0.39$ | $-0.11\pm 1.33$ $-0.09\pm 1.26$ |

Finally, the global RMS of the Zernike coefficients obtained by fitting the irregular part of the corneal topography is summarized in Table 5 . Average RMS values for the different groups are given for the pre- and post-LASIK monozone models, as well as for the optical zone alone. LASIK surgery seems to increase the RMS values, which means a higher surface irregularity, except for the PlanoScan group, where initial values were high already. Interestingly, postsurgical variability decreases, both within and between groups, which suggests that laser ablation yields more homogeneous patterns of surface irregularities. On average, this surgery does not seem to increase surface irregularities. Among groups, Allegretto seems to produce a higher relative increase in Zernike coefficients, but the final values seem to be within normal values. The Zernike RMS values for the post-LASIK OZ are much lower because this optical zone is more regular, but also because the normalization diameter used to compute the coefficients was different: $6\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ for the optical zone and $9\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ for the monozone models. In OZ we can appreciate equivalent, but more marked, differences among groups.

## Table 5

RMS of the Zernike coefficients (in micrometers) resulting after model fit to the corneal topographies (monozone for Pre-LASIK and Post-LASIK, and multizone for OZ).

Pre-LASIK | Post-LASIK | OZ | |
---|---|---|---|

PlanoScan | $9.34\pm 5.60$ | $8.77\pm 3.12$ | $1.11\pm 0.41$ |

Zyoptix | $6.40\pm 3.51$ | $7.65\pm 2.69$ | $0.743\pm 0.19$ |

Allegretto | $6.15\pm 2.91$ | $9.06\pm 1.97$ | $1.40\pm 0.29$ |

Average | $7.64\pm 4.61$ | $8.37\pm 2.67$ | $1.03\pm 0.40$ |

## 4.

## Discussion

So far, we have presented a multizone model of corneal topography for postsurgical eyes. Here it is essential to have the correct number of zones as well as the area covered by each zone. The number of zones was chosen a priori according to the type of surgery. In the particular case of myopic LASIK, the model considers three: OZ, TZ, and P. However, the size of the zones (masks) resulting from the segmentation algorithm differ from the programmed ones ( $5.50\text{-}\mathrm{mm}$ vs. $6.57\text{-}\mathrm{mm}$ average values). The programmed diameter of the OZ is a nominal value, whereas the resulting effective diameter is estimated through a segmentation in which points are assigned to a given zone according to the homogeneity of the descriptors. To determine the source of this bias, we conducted a series of control tests, changing parameters and descriptors in the segmentation algorithm. We found that the segmentation was highly robust against changes in parameters. However, we obtained significantly worse results when removing one or two descriptors. We tried different segmentations using only two parameters (curvature-fit error or curvature-distance), but the results were less robust and reliable, and the resulting zones often showed several disconnected areas.

Mainly, descriptor $D$ (distance to the center) showed a significant effect in decreasing the effective optical zone. In subject AFD, for example, the diameter of the effective optical zone was $5.46\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ versus a $6.5\text{-}\mathrm{mm}$ programmed nominal value. When we repeated the segmentation, removing descriptor $D$ (that is, using only two descriptors), the result was $5.93\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ . Therefore, descriptor $D$ appears to introduce a bias of about $0.47\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ on the estimation of the effective optical zone. Nonetheless, there is a remaining offset of $0.57\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , which cannot be accounted for by segmentation biases. We conclude that from the $1\text{-}\mathrm{mm}$ average difference, about 50% is due to segmentation bias, whereas there is another 50% contribution that is an effective overlapping by the transition zone in terms of lack of homogeneity of descriptors (Gauss curvature and fit error). This sort of uncertainty in the limits of effective zones reveals that there is a statistical uncertainty, but also that it is difficult to establish unique criteria and metrics to determine the effective boundaries. In fact, a major reason to include a transition zone in LASIK treatments is to avoid discontinuities, and in practice this means to avoid prints or marks (i.e., no real physical boundaries). Conversely, nominal programmed zones may differ from effective ones due to uncertainties associated to alignment, eye movements, laser spot size, etc. Therefore, the definition of effective zones is an open question that deserves further study. Descriptor $D$ , proposed here, seems to cause an underestimation of the diameter of the estimated optical zone, but it guarantees a robust and efficient segmentation, which in turn provides good reconstructions of post-LASIK topographies. Nevertheless, the proposed segmentation method is flexible and robust enough so that the definition of descriptors and metrics could be further optimized to avoid (or to minimize) biases, which will be the subject of future work.

The comparative study revealed considerable differences among the three types of treatment. Initial differences among the groups of patients, such as somewhat higher levels of myopic correction in the Allegretto group, or some more regular corneas in the Zyoptix group, may explain some, but not all, differences. Geometrical parameters such as conic constants are known to have a decisive impact on spherical aberration. Figure 10 compares the average wavefront aberration for the pre-LASIK and for the three post-LASIK groups. These aberrations were computed by ray tracing on the topography model for a $6\text{-}\mathrm{mm}$ pupil diameter using Matlab code. Error bars represent standard deviations among subjects. Three main conclusions can be extracted from this figure:

1. There are significant differences in the outcomes of the different treatments, but, surprisingly, differences do not seem related to customization in this sample of patients. Zyoptix (custom) and PlanoScan (standard) yield totally equivalent outcomes with no statistically significant differences. However, corneas treated with Allegretto (custom) show significantly higher RMS wavefront error. This is always the case either if we consider total wavefront (averages: Allegretto $2.27\pm 0.40\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ , Zyoptix $1.60\pm 0.72\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ , and PlanoScan $1.62\pm 0.56\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ ), HOA (Allegretto $1.71\pm 0.34\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ , Zyoptix $1.01\pm 0.48\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ , and PlanoScan $1.02\pm 0.35\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ ), third-order aberrations (coma and trefoil), or fourth-order spherical aberration (Allegretto $1.31\pm 0.20\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ , Zyoptix $0.81\pm 0.39\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ , and PlanoScan $0.74\pm 0.50\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ ). Only second-order aberrations (defocus and astigmatism) are similar for the three treatments.

2. HOA increase in all groups of patients. For Zyoptix and PlanoScan, HOA increase by a factor of about two on average, whereas for the Allegretto group, HOA increase by a factor of three approximately. Results in Fig. 10 suggest that most of this increase is due to spherical aberration, but third-order aberrations (coma and trefoil) also increase. In particular, the increase in third-order aberrations is almost comparable to that of spherical aberration for the Allegretto group. The relative contribution of spherical aberration to the higher-order aberrations is important for the three groups. Before surgery, spherical aberration represents about 50% of HOA, whereas after surgery 80% of HOA is due to spherical aberration. Customization of the treatments does not seem to help reduce spherical aberration in this particular sample of patients.

3. The improvement in second-order aberrations seems modest apparently. However, the defocus coefficient is meaningless for the cornea alone, whereas corneal astigmatism improves after LASIK in all cases. The Zernike coefficient, ${Z}_{5}$ (horizontal-vertical astigmatism), decreases significantly for all of the operated corneas, from an initial average value of $-1.12\pm 0.89\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}-0.41\pm 1.37\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ , $-0.90\pm 0.79\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ , and $-0.54\pm 0.70\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ for Allegretto, Zyoptix, and PlanoScan, respectively.

In summary, HOA increase considerably after surgery.^{28} PlanoScan (standard) and Zyoptix (custom) show totally equivalent outcomes, consistent with previous findings,^{29} namely that HOA increase by a factor of about two after standard surgery. Custom surgery did not seem to improve these results, at least in this particular group of patients. The increment of SA is consistent with the change in the conic constants, which pass from negative to positive values after surgery. These results suggest that the optical quality of the cornea deteriorates with surgery in both custom and standard treatments. In this sense, PlanoScan and Zyoptix provided totally equivalent outcomes, whereas HOA increased more in the Allegretto group, where the levels of refractive correction were higher.

In conclusion, the multizone model permitted us to extract and analyze more and more precise information. It is true that much of the analysis of the optical zone can be done on a monozone model by considering the central 6 or $6.5\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$ , but the multizone model adds further valuable information. First of all, and this was the primary goal of the method, by considering different zones we can keep a high fidelity in the topography model after surgery. Results in Fig. 7 show that the RMS fit error (fidelity metric) obtained with the multizone model after surgery is equivalent to that of the pre-LASIK model. In addition, the RMS error seems slightly lower and more homogeneous both between and within (lower error bars) groups. On the other hand, this type of model allows us to analyze the morphology and alignment of effective zones, defined in terms of homogeneity of the descriptors used rather than nominal values. The set of three descriptors used in the present implementation provided reasonable results, but these are not unique and perhaps better ones could be used in future work. In the comparative study, we focused mainly on the optical zone, but an independent analysis of all the effective zones might help gain further insight on the effects of laser ablation.

## Acknowledgments

This work was supported by the Ophthalmologic Clinic RealVisión and by the Comisión Interministerial de Ciencia y Tecnología (Spain), under Grant FIS2005-05020-C03-01.

## References

**,” Refract. Corneal Surg., 5 (6), 368 –371 (1989). 1042-962X Google Scholar**

*Methods of analysis of corneal topography***,” J. Opt. Soc. Am. A, 23 (2), 219 –232 (2006). https://doi.org/10.1364/JOSAA.23.000219 0740-3232 Google Scholar**

*Optics of the average normal cornea from general and canonical representations of its surface topography***,” Ophthalmic Physiol. Opt., 13 68 –72 (1993). 0275-5408 Google Scholar**

*Mathematical models of the general corneal surface***,” Optom. Vision Sci., 72 821 –827 (1995). https://doi.org/10.1097/00006324-199511000-00008 1040-5488 Google Scholar**

*A spline surface algorithm for reconstruction of corneal topography from a videokeratography reflection pattern***,” J. Opt. Soc. Am. A, 12 2105 –2113 (1995). 0740-3232 Google Scholar**

*Representation of videokeratoscopic height data with Zernike polynomials***,” IEEE Trans. Biomed. Eng., 48 87 –95 (2001). https://doi.org/10.1109/10.900255 0018-9294 Google Scholar**

*Optimal modelling of corneal surfaces with Zernike polynomials***,” Invest. Ophthalmol. Visual Sci., 44 4676 –4681 (2003). https://doi.org/10.1167/iovs.03-0190 0146-0404 Google Scholar**

*Zernike polynomial fitting fails to represent all visually significant corneal aberrations***,” Optom. Vision Sci., 82 1038 –46 (2005). 1040-5488 Google Scholar**

*Automated decision tree classification of corneal shape***,” Invest. Ophthalmol. Visual Sci., 35 2749 –2757 (1994). 0146-0404 Google Scholar**

*Automated keratoconus screening with corneal topography analysis***,” Jpn. J. Ophthalmol., 50 409 –416 (2006). 0021-5155 Google Scholar**

*Automated keratoconus detection using height data of anterior and posterior corneal surfaces***,” J. Refract. Surg., 22 539 –545 (2006). 1081-597X Google Scholar**

*Corneal higher order aberrations: A method to grade keratoconus***,” J. Cataract Refractive Surg., 25 1069 –1079 (1999). 0886-3350 Google Scholar**

*Corneal topography classification in myopic eyes based on axial, instantaneous, refractive, and profile difference maps***,” Appl. Opt., 44 4528 –4532 (2005). https://doi.org/10.1364/AO.44.004528 0003-6935 Google Scholar**

*Differences between real and predicted corneal shapes after aspherical corneal ablation***,” Ophthalmology, 110 1926 –1930 (2003). 0161-6420 Google Scholar**

*Difference map or single elevation map in the evaluation of corneal forward shift after LASIK***,” J. Cataract Refractive Surg., 31 2350 –2355 (2005). 0886-3350 Google Scholar**

*Goodness-of-prediction of Zernike polynomial fitting to corneal surfaces***,” IEEE Trans. Biomed. Eng., 51 2203 –2206 (2004). 0018-9294 Google Scholar**

*A refined bootstrap method for estimating the Zernike polynomial model order for corneal surfaces***,” Modern Differential Geometry of Curves and Surfaces with Mathematica, 373481 –380500 2nd ed.CRC Press, Boca Raton, FL (1997). Google Scholar**

*The Gaussian and mean curvatures and surfaces of constant Gaussian***,” J. Cataract Refractive Surg., 31 1607 –1613 (2005). 0886-3350 Google Scholar**

*Reliability of Orbscan II topography measurements in relation to refractive status***,” Ophthalmology, 113 177 –183 (2006). 0161-6420 Google Scholar**

*Evaluation of Orbscan II corneal topography in individuals with myopia***,” J. Cataract Refractive Surg., 31 379 –384 (2005). 0886-3350 Google Scholar**

*Functional optical zone after myopic LASIK as a function of ablation diameter***,” Invest. Ophthalmol. Visual Sci., 48 1053 –1060 (2007). 0146-0404 Google Scholar**

*Functional optical zone of the cornea***,” J. Cataract Refractive Surg., 28 948 –953 (2002). 0886-3350 Google Scholar B. S. Boxer Wachler, V. N. Huynh, A. F. El-Shiaty, and D. Goldberg, J. Cataract Refractive Surg., 28 942 –947 (2002). https://doi.org/10.1016/S0886-3350(02)01324-X 0886-3350 Google Scholar**

*Evaluation of corneal functional optical zone after laser**in situ*keratomileusis**,” J. Cataract Refractive Surg., 29 762 –768 (2003). https://doi.org/10.1016/S0886-3350(02)01895-3 0886-3350 Google Scholar**

*Changes in corneal asphericity after laser**in situ*keratomileusis**,” J. Cataract Refractive Surg., 29 2096 –2104 (2003). https://doi.org/10.1016/j.jcrs.2003.09.008 0886-3350 Google Scholar**

*Spherical aberration after laser in situ keratomileusis and photorefractive keratectomy—clinical results and theoretical models of etiology***,” J. Cataract Refractive Surg., 31 127 –135 (2005). 0886-3350 Google Scholar**

*Causes of spherical aberration induced by laser refractive surgery***,” J. Opt. Soc. Am. A, 23 544 –549 (2006). 0740-3232 Google Scholar**

*Spherical aberration of the anterior and posterior surfaces of the human cornea***,” J. Cataract Refractive Surg., 27 362 –369 (2001). https://doi.org/10.1016/S0886-3350(00)00806-3 0886-3350 Google Scholar**

*Increased higher-order optical aberrations after laser refractive surgery: a problem of subclinical decentration***,” Invest. Ophthalmol. Visual Sci., 42 1396 –1403 (2001). 0146-0404 Google Scholar**

*Ocular aberrations before and after myopic corneal refractive surgery: LASIK-induced changes measured with laser ray tracing*