Measurement, modeling and visualization of the shape of the cornea
F.M. Vos et al.

[Abstract], [Introduction], [Image Processing], [Surface Reconstruction],
[Conclusions and Future Work], [References]

Section 3: Surface Reconstruction

A: Introduction

The goal of corneal topography algorithms is to determine the shape of the corneal surface. In the past a number of algorithms have been proposed to convert the input information to a (mathematical) description of the cornea (klein94, andersen93, halstead95). At first all techniques assumed axial symmetry resulting in a reconstruction of the cornea in a finite number of cross cuts. Characteristic for most of these latter methods is that given a position i on the cornea (corresponding to a ring transition) the shape of the corneal profile between the current position and a next (i+1) is reconstructed. As a consequence errors progress from one ring to the next. In andersen93 a new method was presented that estimated the shape of the cornea with an iterative procedure using a global ellipse model.

A big step forward has been made by halstead95 who provided an algorithm modelling the cornea as a B-spline surface. Halstead also models the shape of the eye by iterative approximation. However, as he developed his method using a (classical) Placido stimulator, assumptions had to be made as to the unique correspondence of positions on the CCD and on the stimulator.

Our approach allows for an integral surface reconstruction of the corneal shape withour any prior assumptions.

B: General Description

Our objective is to model the shape of the eye as a parameterized surface: z = T(theta, x, y). Suppose a point s on the stimulus is reflected of the cornea in a point c on the CCD-chip as shown in this diagram
Given an instance of our model-function, it is possible to calculate a point s' on the stimulator that would be registered in c when mirrored to the surface described by T. This can be done by simple backwards ray-tracing. Our objective is to minimize || s' - s || ² for the crossings in our checkerboard pattern, coming down to nonlinear least squares parameter-estimation of the parameters (theta) of the model-function. To this aim standard parameter-estimation procedures are available.

C: Spline Theory

For a first approximation of the shape of a normal eye an ellipsoid could be taken as a model. Possible aberrations of this limited shape could exist by way of very local bumps or more global bulging of the surface. To model these aberrations other types of surfaces are required with more degrees of freedom.
One model that fullfills this requirement is introduced in halstead95. The so-called tensor product B-spline with uniformous knot spacing employed is build out of small "pieces" of surface: patches. Each patch is described by the coefficients of a polynomial of preset degree. The degree of the polynomial is determined for the largest part by the information it should represent.
As physicians are mainly interested in curvature information of the surface second derivative), the surface should at least consists of patches of degree 2 (in the latter case the curvature is represented by a constant value). We have opted for splines of degree 5 with which the curvature can be represented by polynomials of degree 3.

The parameters that determine the shape of the spline are the so-called control-points. A spline surface of m by n patches is defined by (m + 5) x (n + 5) control-points. The surface is constructed as the sum of scaled basis functions. The control points (c_ij) determine the scale of the basis functions:

S(x, y) = Sum _i { Sum _j { c _ij b _i (x) b _j (y) } }

where i and j run from 1 to m+5 resp. n+5. With the type of splines with "uniformous knot spacing" each basis function is the translated copy of the standard basis function b(u) given in this figure. The spline basis function is non zero in 6 unit intervals, within which it is defined (in our case) as a fifth degree polynomial. The scaled basis functions are defined as: b_i(x) = b(x-i-1) and b_j(y) = b(y-j-1). The following figure gives two examples of the 1-dimensional analogon of this type of spline. From the previous it can be concluded that in principal S(x,y) can be evaluated in every arbitrary position (x,y). However, only in a restricted area of the domain S evaluates to non zero (depending on the number of patches).
A patch is defined as unit square within which the surface is defined by the sum of exactly 36 basis functions. Commonly , the domain of S is restricted to the area in which non less than 36 basis functions evaluate to a value unequal to zero.
A slight variation on this first type of B-splines is characterized by nonuniformously distributed basic functions that do not have identical shapes (in computer graphics terms: bivariate tensor product B-splines with nonuniformous knot spacing). This figure shows examples of this type of B-spline.

D: Results

As an evaluation of the measurement technique, the algorithm was tested using six objects of known shape: three spheres (with R=6.5, 7.0 respectively 8.0 mm) and three rotationally symmetrical ellipsoids of which the shape was measured by mechanical techniques to obtain a gold standard.
A B-spline build from one patch was used for a model. The next table gives the RMS-error of the measurements of the spheres as well as the average curvatures of the respective models and the corresponding standard error
The RMS error is defined as the root mean square between the calculated and the ideal surface, corrected for a systematic positional bias. From the table it can be concluded that within micron precision the surfaces are recovered.

**Results for spherical test objects**
R object	RMS-error (mm)	Average Mean Radius + std. dev.(mm)
6.5	7.06 10^-4	6.05 (0.03)
7.0	1.02 10^-3	7.03 (0.04)
8.0	8.92 10^-4	7.99 (0.03)

For the three rotationally symmetrical ellipsoids the results are as follows

**Results for ellipsoid test objects**
a/b	RMS-error (mm)
0.99	1.12 10^-3
0.93	1.16 10^-3
0.87	1.13 10^-3

For these symmetrical objects a numerical measure can readily be defined to determine the correctness of the approximation. However with more complicated shapes - as will be the case in pathological interesting cases - such a measure is not readily available. Then visual feedback between the user and the analysis software is of paramount importance.
Again the visualization environment DX is used to provide this feedback. In this paper we will show two relevant cases.

The first case has to do with the effect of the number of patches used for the surface approximation. For test purposes we simulated the response of a spherical object (radius 8mm) on which a small distortion of 2.2 mm diameter and height 0.02 mm was placed (using PovRay) and analyzed the thus defined image with the methodology outlined above.
Four different approximations are used: splines of 1x1, 2x2, 4x4 and 8x8 patches. Rather than displaying the actual surface we have visualized the difference between the actual surface and a perfect sphere.
These visualizations clearly show that as the number of patches increases a better approximation of the surface is obtained or rather that a high number of patches is necessary for irregularly shaped objects. The drawback is that the computational time increases exponentially with the number of patches which makes this kind of approximation highly non-interactive.

Another way of presenting the data is differentiating the calculated surface twice, as a first approximation of a calculation of the curvature and then determining iso-surfaces. Applying this to real data yields the following result . The program is given in this figure.

The use of Data Explorer is relevant at two levels: one is that of platform for the development of software tools. The previous example shows that it is easy to prototype complex situations. Espcially the routine 'compute' is handy in such situations.
For more advanced problems the possibility to integrate via a well defined program interface locally developed C routines with existing code adds to the usability of the product.