![MathJax Logo](/templates/jsp/_style2/_tandf/pb2/images/math-jax.gif)
Abstract
In this paper, we apply an iterative regularization technique based on Newton’s method to the ill-posed problem of parameter identification in a differential equation describing corneal surface. This equation is based on our previously derived model of corneal topography. We prove the convergence of the proposed method along with some stability estimates. When compared with real corneal data, this method gives reasonable results with error of the same order of magnitude as is introduced by measuring equipment in data acquisition. Moreover, the iteration converges very quickly to the parameters associated with the intra-ocular pressure and the elasticity of the cornea.
1 Introduction and model description
In present times, ophthalmology is one of the most rapidly developing fields in medical sciences. There is a need for a better treatment of various eye diseases since increasingly more people develop sight problems. This is mostly due to the ‘civilization’ issues present in our lives, i.e. excessive reading and computer use. Since sight is almost by anyone considered as the most important sense, keeping healthy eyes is a crucial matter. One of the most important parts of human eye is the cornea. It is a transparent, shell-like structure that constitutes a large portion of the frontal part of the eye. Together with sclera, it constitutes the outer part of the eyeball. Cornea has a layered structure with its thickness of about mm. It is responsible for about two-thirds of refractive power and, partly because of this, so much attention is being put in studying its nature.[Citation1] Moreover, cornea is very durable and acts as a frontline defence against the exterior conditions. For more on biology of cornea see [Citation2]. The beginnings of mathematical treatment of the cornea dates as back as to Helmholtz [Citation3]. He was an inventor of ophthalmoscope, investigated optics of seeing and colour perception. For a modern mathematical consideration concerning a broad range of eye-associated issues, the reader is invited to see [Citation4, Citation5]. There, fluid-mechanical models of different flows and their effect on sight disorders are proposed. The issues of tear-film dynamics and flows are considered for example in [Citation6–Citation8]. Particularly, Authors develop a model of so-called break-up time (BUT) which is the first moment after a blink of appearance of a dry patch on the cornea. The cited papers show that mathematics which is needed in modelling eye phenomena is far from trivial.
The development of medical equipment and new means of examination gave the need for creating mathematical descriptions and models of various phenomena associated with sight. Particularly, the shape and topography of the cornea are very important factors used in diagnosing many sight disorders. The most common corneal shape models that are presently in use are based on surfaces of revolution generated by conic sections (see [Citation9, Citation10]). Unfortunately these models, although very useful and accurate, lack the physical motivation what makes them conceptually not fully justified. In that case, an usual attempt on modelling cornea is made by creating the best-fit ellipsoid to the surface of the cornea. The calculated shape parameters are then used in further reasoning (see e.g. [Citation11]). Another group of corneal topography models consists of those based on shell theory and/or finite-element method (e.g. [Citation12, Citation13]). These models are sometimes too complex to be readily computed and used in real-world applications. Nevertheless, they can serve many purposes such as reconstruction of the corneal surface by the eye measurements. A large number of scientists use the Zernike orthogonal polynomials to describe the shape of the cornea (see e.g. [Citation14, Citation15]). These polynomials are orthogonal with respect to a certain scalar product and can be used to describe aberrations in (organic) lens. The advantage of these polynomials is that they are indexed by parameters which have a straightforward physical meaning (see [Citation16, Citation17]). Apart from Zernike, Bessel functions have also been used for similar aim.[Citation18] Also, adaptive techniques of approximating corneal surfaces are also being used.[Citation19]
In [Citation20] we have proposed a new mathematical model of a human cornea. This model was further improved in [Citation21]. We have derived the following PDE that governs the elevation of the cornea over its domain
1
1 where
and
are nonnegative, dimensionless constants related to the tension
, intra-ocular pressure
and elasticity constant
through
2
2 Here,
and
are dimensionless quantities scaled with respect to the typical corneal radius
(which is on the average of about
mm). The boundary condition says that the cornea has a rim at
. This PDE has been derived by describing the corneal surface as a thin membrane with constant tension
, elastic force proportional to the elevation
and a force caused by the intra-ocular pressure
acting normally to the surface. Assumption that
and
are constant (hence the pressure is also constant throughout the whole cornea) suggests that we should look for an axisymmetric solution. This setting can be thought as a zero-order perturbation approximation to the real situation, where cornea is asymmetric and aspherical. This was partially justified and investigated in [Citation22], where the inverse problem associated with variable
was elaborated and shown that corrections to the axisymmetric topography are of higher order. Another simplifying assumption we make is that the pressure–force acts along the
-direction and not in the direction of the normal vector. From the derivation in [Citation21], it can be seen that this clears off the square root on right-hand side of (Equation1
1
1 ) (square root is a cosine of the angle subtended between
-axis and normal vector). This is not a very severe assumption, since we can assume that the corneal radius of curvature is small. This can be tentatively justified by a crude estimate: typical elevation of the cornea is about
mm and its radius
mm. This gives
, which implies that
3
3 which is close to
.
With , we get an ODE for the corneal cross-section:
4
4 with boundary conditions
5
5 The condition
is obvious since we want cornea to have a rim at
and
gives us a smooth at origin, nonsingular surface of revolution. Equation (Equation4
4
4 ) can be readily solved yielding
6
6 where
is a Modified Bessel Function of the First Kind. We note that when fitting with the real corneal data we have to consider both interior and exterior corneal surfaces. This is due to the fact that cornea has a thickness. Also, we remark that in the future work, thickness will be incorporated into our differential equation model.
An easy calculation can give us the estimate on the error that is made when neglecting the square root. In [Citation21] it was shown that the solution to the problem (Equation44
4 )–(Equation5
5
5 ) can be stated in terms of the integral equation
7
7 while (Equation6
6
6 ) can be rewritten as:
8
8 where
is a Green’s function for our problem. Taking the difference we can estimate
9
9 Thus making the assumption of pressure acting vertically introduces a maximal error of
. This is an a-priori theoretical estimate. We can make a more realistic calculation by utilizing the typical corneal sizes and (Equation3
3
3 )
10
10 which would give us
mm, where we took the typical corneal deflection to be equal to
mm. For the accuracy of the first approximation we can thus safely neglect the square root in (Equation4
4
4 ).
In this paper, we are mainly concerned with determining constants and
in (Equation6
6
6 ) when the solution
is known (usually with a noise). This is one example of so-called inverse problems (as opposed to the Direct Problems, here: solving the ODE (Equation4
4
4 )). In applications one usually wants to make an identification of various parameters in differential equations (see [Citation23]). In our case, knowing
and
will give us an approximation to the corneal topography. A problem of finding (possibly nonconstant)
when
is known was solved in [Citation22].
Inverse problems are usually ill-posed, that is they can have no solutions, have a nonunique solution or be unstable with respect to the noise level in input data.[Citation24] In our case we cannot expect that the cornea is mathematically described precisely by (Equation66
6 ) hence generally our inverse problem will have no solution. The simplest way of determining unknowns
and
in (Equation6
6
6 ) is by using maximum deflection
and central radius of curvature
as was done in [Citation21]. However, this method is not reliable because
depends on a specific reference level and could be determined ambiguously. Moreover, it can give serious errors because of instability.
Due to the nonuniqueness there is a need of defining some generalized solution to the inverse problem. This is usually done by looking for the unique least-squares solution.[Citation25] Moreover, we often do not know the exact function . When doing measurements some noise
would always be present, thus we have only
at our disposal that is,
11
11 It is natural to expect a stable dependence of parameters
and
on the noise level
. Small errors in measurement should give small errors in parameters. In parameter identification problems usually this is not the case and we have serious instabilities in the obtained solutions. Such parameters are not of much use in applications.
We shall put our inverse problem in a slightly more abstract setting of space and propose a method of determining the unique solution to the inverse problem of the type we are concerned of. To overcome these problems, many regularization methods have been developed (see e.g. [Citation26]). We present a variation of one of them which gives us a stable and unique solution to our parameter identification problem.
2 Analytical results
In this section, we propose an iterative regularization method and prove some stability results. First, we introduce the following family of linear operators defined by12
12 where
is a fixed function with possibly nonlinear dependence on
. Here we assume that
(
compact) and
, hence
13
13 Moreover, we assume that the family of operators
depends smoothly (up to the second order) on
.
We treat the real line as a Hilbert space with multiplication as a scalar product. The adjoint
of
is defined by
14
14 for all
. Thus, by (Equation12
12
12 ), we have
. It is straightforward to observe that
15
15 where appropriate types of norms are understood.
By the inverse problem we will understand the determination of and
from the equation
16
16 where
is a known function (in our case
will describe the corneal geometry and thus
must satisfy the boundary condition (Equation5
5
5 )). It is not expected that the cornea will be exactly described by the function
that is,
does not necessarily have to be in the range of
. Moreover,
is always measured with noise
which also has to be taken into account. To make precise what we understand by an inverse problem, we must first make a distinction between the well-posed and ill-posed problems (in the Hadamard sense, see e.g. [Citation27]).
Definition 2.1
A mathematical problem (for example a differential or algebraic equation) is called a well-posed problem if it has all of the following properties
it has a solution,
it has only one solution,
it is continuous with respect to the perturbations in the initial data (stability).
If a problem is not well-posed it is called ill-posed. Inverse problems are usually ill-posed thus are more difficult to analyse and use in applications than the well-posed ones.
The Equation (Equation1616
16 ) is linear in
and generally nonlinear in
. Assume for now that
is fixed and consider only finding
as a solution to the linear equation. We look for a generalized least-squares solution, that is
is such solution when
17
17 Of course, when
is an exact solution to (Equation16
16
16 ), the left-hand side of (Equation17
17
17 ) is zero, thus it is also a least-squares solution. It can be shown (see e.g. [Citation25]) that the sufficient and necessary condition for
to be the unique least-squares solution to (Equation16
16
16 ) is
18
18 From (Equation18
18
18 ) using (Equation12
12
12 ) and (), we obtain
19
19 It is evident that this least-squares solution depends on parameter
and
is uniquely determined when
is known. Now, the solution to the full nonlinear problem (Equation16
16
16 ) is reduced to finding only one parameter, namely
. To simplify notation let us define a nonlinear operator
20
20 and consider finding least-squares solution
to the problem
. Due to the compactness of
, the existence of this solution is provided. We will proceed by a method known as Output Least Squares (see [Citation23]). This iterative method takes from Newton’s tangent method of finding roots of nonlinear equations. The difference is that we look for a generalized least-squares solution and not the exact one.
Fix an initial guess of unknown parameter . We will try to find an increase
such that
. We can write
21
21 where
is a Fr
chet derivative, which is a linear operator and remainder
. In our simple case of one-dimensional operator (Equation12
12
12 ), the Fr
chet derivative is just
22
22 It is easy to see that the adjoint
has a form
23
23 with the following equalities in norms
24
24 We look for
as a least-squares solution to (Equation21
21
21 ), thus discarding remainder we want to have
25
25 Now, using (Equation22
22
22 ), () and (Equation24
24
24 ) we get
26
26 Using
as a new initial guess, we can continue this procedure and propose an iterative scheme
27
27 Proposed algorithm can be summarized in the following lines, where all the formulas are consequences of (Equation19
19
19 ), (Equation20
20
20 ), (Equation12
12
12 ), () and (Equation27
27
27 ).
(1) | Let | ||||
(2) | Chose a starting value | ||||
(3) | Compute | ||||
(4) | Let | ||||
(5) | Repeat (3)–(4) until the desired accuracy |
Theorem 2.2
Assume that28
28 then there exists a neighbourhood of
in which the iterative scheme (Equation27
27
27 ) is convergent to the least-squares solution of
.
Remark 1
By , we mean a bilinear functional
. It is straightforward to see that
. Thus we write
. We also will make an abbreviation by writing
instead of
.
Remark 2
The assumption (Equation2828
28 ) is not very severe because we expect that the residue
would usually be small (for the method to be useful). Of course in applications we do not know
, but sometimes we can at least try to guess its approximate location. Thus, due to minimality of
, we can impose a stronger assumption
29
29 which implies (Equation28
28
28 ).
Proof
Denote an error of approximation by . We will show that
. We have
30
30 Observe that
31
31 hence application of
gives
32
32 where
is between
and
. We can rewrite (Equation30
30
30 ) as
33
33 We are looking for the least-squares solution
, thus the normal equation
must be fulfilled (see [Citation25]). This is the particular difference between the classical proof of Newton’s tangent method, where the exactness of the solution is utilized. Therefore, writing
34
34 we finally arrive at
35
35 which gives
36
36 From (Equation28
28
28 ) we have
. Thus, when
is sufficiently small, that is, if we start our iteration suitably close to
, the assumption (Equation28
28
28 ) implies that
37
37 where
. This means that
which is what we had to show.
Now, we have a method of determining and
as a least-squares solutions to
. Next, we will show some stability results. When we know only noisy data
, the iteration is convergent to some
which usually is different from
. This produces an error bound which grows with
. The question arises in what moment should we stop iteration (Equation27
27
27 ) to have the smallest bound possible. The next Theorem gives an a-priori stopping rule.
Theorem 2.3
Assume we have a noisy data , then for solutions
to
we have the following estimates
38
38 where
is some
-dependent constant and
is an optimal stopping moment chosen to satisfy
39
39
Proof
Let and
be sequences (Equation27
27
27 ) convergent to
and
, respectively. We will use the following fundamental estimate
40
40 where the error decomposes into two terms: one corresponding to the noise and the other connected with the convergence rate of the iterative method. We can take
and estimate the first term as follows:
41
41 Using the fact that
and
are Lipschitz continuous (since they have bounded
-derivatives) with respect to
and that
is a compact set, we have
42
42 where we used the fact that
is bounded with respect to
(due to compactness of
) and put it along other estimates under
. We can further write
43
43 where
is greatest constant of the following
44
44 Using (Equation43
43
43 ) it is easy to see that
and
. By induction we can finally prove that
45
45 The second term in (Equation40
40
40 ) is by (Equation37
37
37 ) dominated by
, thus we can write
46
46 We see that with increasing
one error bound grows without a limit and the other goes to zero. Thus, in order to provide its convergence with
we must find a minimizing
as a function of
. Simple calculations give this minimizing argument equal to
47
47 Plugging this expression into the estimate (Equation46
46
46 ) we obtain
48
48 By some elementary manipulations along with the identity
we arrive at
49
49 whereby
and
we denoted the
-independent constants. We see that the exponent of
is between
and
. For the estimate of
denote
and by (Equation19
19
19 ) write
50
50 which gives us the desired estimate for the second parameter and concludes the proof.
3 Numerical simulations and application to the corneal geometry
In this section, we apply our previous abstract reasoning to the problem of finding and
in the Equation (Equation4
4
4 ) describing the corneal geometry. We will find a fit with the two data-sets consisting of
points representing the interior and exterior surfaces of cornea. Unfortunately, we are in possession of this data only. A more thorough verification would only be possible provided a larger set of data.
Instead of (Equation66
6 ) we will use the following function
51
51 which is
in (Equation6
6
6 ). This is done because of the form in which data are acquired from the measuring equipment. In the case of our measurement, the noise level is
m. Of course, we have to return to the original dimensional model, but here we suppress the details.
Our previous work [Citation21] gives us a hint where to look for the unknown parameters and
, where we tentatively found that
,
,
and
. We expect that the real values of parameters will be of order of unity, thus we restrict our search of
for example to the set
. Numerical calculations show that both
and
are decreasing functions of
.
We note that the needed assumption of Theorem 2.2 is satisfied. Indeed, we have52
52 for exterior and interior surfaces, respectively. These quantities are less than
so the iteration (Equation27
27
27 ) is convergent. For interior and exterior surfaces we have, respectively,
,
and
,
. The number of iterations needed for the convergence for different starting points is presented in Table . We can see that even for big differences between the starting point and the limit, only a few iterations are needed. When doing the same calculations with other iterative regularization method, namely Landweber iteration (see e.g. [Citation26]), achieving the same precision took almost
and
iterations, respectively! However, Landweber’s iteration is known to converge arbitrarily slow (see [Citation28]). Our results were verified by using some well-known algorithms for finding the least-squares solution. We have used Newton, Quasi-Newton, Levenberg-Marquardt and Conjugate Gradient methods available in the scientific environment MATHEMATICA. In any case we have obtained the same values of
and
.
Table 1. Number of iterations needed for the convergence with different starting points for interior and exterior surfaces of the cornea.
Table 2. Results of the numerical experiment for finding true values of parameters and
. The model was constructed from (Equation51
51
51 ) with additive noise on level
and
,
.
By Theorem 2.3 and by using the known value of , we crudely estimate that the optimal number of iterations should be about
. Also, by (Equation38
38
38 ) the difference between the computed
,
and the exact values of the unknown parameters
and
can be estimated to be of order of
. The absolute fitting error was calculated from the formula
53
53 which is depicted in Figure . The plots show the contour plots of the magnitude of the absolute error
along the region of cornea. The means of these errors taken along the surface of the cornea are equal to
mm for the interior surface and
mm for the exterior. These errors are within the systematic measurement uncertainty, which is
m. Therefore we do not introduce any additional errors which are of lower order than the data. The dimensional radius of the cornea is
mm and the reference surface of the measurement is drawn in light colour on the outside of the circular rim of the cornea (that is for
mm).
Figure 1. Contur plots of the absolute errors in fitting defined in (Equation53
53
53 ). On the top: exterior surface, on the bottom: interior surface of cornea. Horizontal and vertical axes, shown in millimetres, depict
and
coordinates of the point on the corneal surface. The light colour outside the cornea represents the reference surface.
![Figure 1. Contur plots of the absolute errors in fitting E(x,y) defined in (Equation5353 E(x,y):=h(x,y)-F(an∗†)(x,y),53 ). On the top: exterior surface, on the bottom: interior surface of cornea. Horizontal and vertical axes, shown in millimetres, depict x and y coordinates of the point on the corneal surface. The light colour outside the cornea represents the reference surface.](/cms/asset/4a35d063-e639-474d-8aeb-4332569e3a2e/gipe_a_922074_f0001_b.gif)
Our algorithm can also be tested for a performance with different noise levels . For this numerical experiment, we produce an artificial corneal surface described by the function (Equation51
51
51 ) with an additive random noise generated from normal
distribution. We choose
and
as the true parameters in our simulation. The results of the experiment are presented in Table . The number of iterations for different starting points chosen from the set
is small, ranging from
to
. It is evident that even for a large (comparing to the maximum corneal deflection
mm) noise level, the method shows decent stability yielding results that are very close to the exact parameters.
4 Conclusion and discussion
In this paper, we have applied a simple iterative regularization method to the problem of parameter identification in a differential equation describing corneal surface. At the beginning, we have reduced a double linear–nonlinear problem to a nonlinear one and applied Newton’s iteration to solve it. Sufficient conditions for convergence were proved. We have also found estimates for order of stability with respect to noise in input data.
When fitting with the real data, we saw that Newton’s method needed fewer hundred iterations that Landweber’s. The fitting error is very reasonable and is smaller for exterior than interior surface of the cornea (we have observed similar results earlier in [Citation21]). We have calculated fitting parameters and
which by Theorem 2.3 converge to the real least-squares solutions
and
when
. Thus we have obtained the unique, stable solution to our inverse problem.
In our model, parameters and
are associated with physical quantities (Equation2
2
2 ) such as intra-ocular pressure. Stable method of obtaining them is crucial in the real-world applications since every measurement introduces some noise. Intra-ocular pressure is an important parameter for the diagnostics of many seeing disorders and in almost every eye examination its measurement has to be taken. Since parameter
is associated with the intra-ocular pressure, our model might help to determine it in a noninvasive way, knowing only the corneal topography
. Nevertheless, we hope that our model will provide some new insights into the biomechanical structure of cornea and make a step to determining intra-ocular pressure without any invasive techniques.
Acknowledgements
The authors would like to thank Dr Robert Iskander from Institute of Biomedical Engineering and Instrumentation, Wroclaw University of Technology, Poland and School of Optometry, Queensland University of Technology, Australia for access to the data. This research was supported by the Polish Government funds for science: research project NCN 2012/05/N/ST1/02860.
References
- Mejía-Barbosa Y, Malacara-Hernández D. A review of methods for measuring corneal topography. Optom. Vision Sci. 2001;78:240–253.
- Trattler W, Majmudar P, Luchs JI, Swartz T. Cornea handbook. Slack Incorporated, Thorofare; 2010.
- von Helmholtz H, Southall H, James PC, editors. Helmholtz’s treatise on physiological optics. Dover Publications, New York; 1924.
- Canning CR, Dewynne JN, Fitt AD, Greaney MJ. Fluid flow in the anterior chamber of a human eye. IMA J. Math. Appl. Med. Biol. 2002;19:31–60.
- Ismail Z, Fitt AD. Mathematical modelling of flow in Schlemm’s Canal and its influence on primary open angle glaucoma. International Conference on Science and Technology (ICSTIE): Applications in Industry and Education, 1967–1973; 2008; Universiti Teknologi MARA Penang.
- Braun RJ, Usha R, McFadden GB, Driscoll TA, Cook LP, King-Smith PE. Thin film dynamics on a prolate spheroid with application to the cornea. J. Eng. Math. 2012;73:121–138.
- Braun RJ. Dynamics of the tear film. Annu. Rev. Fluid Mech. 2012;44:267–297.
- Braun RJ, Fitt AD. Modelling drainage of the precorneal tear film after a blink. Math. Med. Biol. 2003;20:1–28.
- Kasprzak H, Iskander DR. Approximating ocular surfaces by generalized conic curves. Ophthal. Physiol. Opt. 2006;26:602–609.
- Rosales MA, Juárez-Aubry M, López-Olazagasti E, Ibarra J, Tepichín E. Anterior corneal profile with variable asphericity. Appl. Opt. 2009;48:6594–6599.
- Kiely PM, Smith G, Carney LG. The mean shape of the human cornea. Optica Acta. 1982;29:1027–1040.
- Anderson K, El-Sheikh A, Newson T. Application of structural analysis to the mechanical behaviour of the cornea. J. R. Soc. Interface. 2004;1:3–15.
- Ahmed E. Finite element modeling of corneal biomechanical behavior. J. Refract. Surg. 2010;26:289–300.
- Iskander DR, Collins MJ, Davis B. Optimal Modeling of Corneal Surfaces by Zernike Polynomials. IEEE Trans. Biomed. Eng. 2001;48:87–95.
- Schneider M. Iskander DR, Collins MJ, Modeling corneal surfaces with rational functions for high-speed videokeratoscopy data compression. IEEE Trans. Biomed. Eng. 2009;56:493–499.
- Schwiegerling J, Greivenkamp JE, Miller JM. Representation of videokeratoscopic height data with Zernike polynomials. J. Opt. Soc. Am. A. 1995;12:2105–2113.
- Schwiegerling J. Cone dimensions in keratoconus using Zernike polynomials. Optom. Vision Sci. 1997;74:963–969.
- Trevino JP, Gómez-Correa JE, Iskander DR, Chávez-Cerda S. Zernike vs. Bessel circular functions in visual optics. Ophthalmic Physiol. Opt. 2013;33:394–402. doi:10.1111/opo.12065.
- Martínez-Finkelshtein A, Ramos-López D, Castro-Luna GM, Alió JL. Adaptive cornea modeling from keratometric data. Invest. Ophthalmol. Visual Sci. 2011;52:4963–4970.
- Okrasiński W, Płociniczak Ł. A nonlinear mathematical model of the corneal shape. Nonlinear Anal.: Real World Appl. 2012;13:1498–1505.
- Okrasiński W, Płociniczak Ł. Bessel function model of corneal topography. Appl. Math. Comput. 2013;223:436–443.
- Płociniczak Ł, Okrasiński W. Regularization of an ill-posed problem in corneal topography. Inverse Prob. Sci. Eng. 2013;21:1090–1097.
- Groetsch CW. Inverse problems in the mathematical sciences. Braunschweig: Vieweg; 1993.
- Hofman B. On ill-posedness and local ill-posedness of operator equations in Hilbert spaces. Fakultät für Mathematik: Technische Universität Chemnitz; 1998.
- Kirsch A. An introduction to the mathematical theory of inverse problems. New York (NY): Springer; 1996.
- Engl HW. Regularization methods for the stable solutions of inverse problems. Surv. Math. Ind. 1993;3:71–143.
- Hadamard J. Le problème de Cauchy et les équation aux dérivée partielle linéaires hyperboliques. Paris: Hermann; 1932.
- Hanke M, Neubauer A, Scherzer O. A convergence analysis of the Landweber iteration for nonlinear ill-posed problems. Numer. Math. 1995;72:21–37.