Measurement 240 (2025) 115613

Contents lists available at ScienceDirect

Measurement
journal homepage: www.elsevier.com/locate/measurement

Partial reconstruction of sparse and incomplete point clouds applied to the

generation of corneal topographic maps of healthy and diseased patients
Alejandro Ballesta a , José S. Velázquez b, Jorge L. Alió del Barrio c,d , Francisco Cavas b,*
a
International School of Doctorate, Technical University of Cartagena, Cartagena, Spain
b
Department of Structures, Construction and Graphical Expression, Technical University of Cartagena, Cartagena, Spain
c
Cornea, Cataract and Refractive Surgery Department, VISSUM, Alicante, Spain
d
Division of Ophthalmology, Department of Pathology and Surgery, Faculty of Medicine, Miguel Hernández University, Alicante, Spain

A R T I C L E I N F O A B S T R A C T

Keywords: The presence of unmeasured areas or areas with missing data in corneal topographies greatly affect the quality of
Corneal model digital corneal models. These are more frequent in posterior than in anterior corneal surface. This work seeks to
Keratoconus complete and improve digital corneal models generated from incomplete point clouds obtained with commercial
Topography
tomographic equipment, applying three interpolation/extrapolation techniques based on population data
Cubic RBF
Biharmonic Splines
(Polynomial, Biharmonic Splines and Cubic Radial Base Functions). 43 personalized corneal models of patients
Polynomic model with different grades of keratoconus were used to test each interpolation method. All interpolation results were
compared using the Mean Average Percentage Error (MAPE) index, observing that it is possible to interpolate/
extrapolate up to the radius 5.2 mm on the anterior face and up to 4.2 mm on the posterior face in general, but
none of the methods stand out above the others for all the situations considered. Consequently, extrapolation
method must be chosen carefully.

1. Introduction in a very short time, obtaining a three-dimensional measurement of the

anterior chamber and its structures. The combination of two mono­
Corneal topography is a non-invasive technique that allows chromatic Scheimpflug cameras with 360-degree rotation and a Placido
measuring the shape and elevation of the corneal surface. This technique disc allows for a complete analysis of the cornea and the anterior
is widely used in clinical practice in the diagnosis and treatment of segment, providing real and reliable data on the elevation of the entire
various eye diseases, such as keratoconus [1,2], as well as in the plan­ anterior segment (referred to as raw data), the tangential and axial
ning of refractive and cataract surgeries [3]. curvature of the anterior and posterior corneal surfaces, and the global
One of the most widely used topographers in clinical practice is the refractive power of the cornea, as well as a biometric estimation of
Sirius topographer (CSO, Italy) because it offers highly precise and various structures and a corneal wavefront map with visual quality
repeatable measurements of corneal curvature, as well as other analysis and corneal pachymetry maps [7].
anatomical measurements of the anterior segment [4]. It is based on the However, this topographer has several drawbacks, such as the time
integration of two technologies: the Placido discs, and the Scheimpflug required to complete the full rotation for data acquisition, which can be
photography principle. The projection of the Placido discs [5] allows affected by extrinsic errors [5], such as tear film stability or errors in
obtaining a specular image by reflection and, through an Arc-Step al­ centering, alignment, or focusing during the measurement acquisition
gorithm, calculating the curvature and the reflection power of the process, or errors due to obstruction of the resulting visual field (eyelash
anterior surface. The projection of a light slit onto the cornea is photo­ shadows) [8]. Another disadvantage of this technology is that it uses
graphed using the Scheimpflug technique [6], whose main characteristic Zernike polynomials for the reconstruction of corneal surfaces (data
is that the camera lens plane is placed not parallel but at an angle to the processed by proprietary internal algorithms), and there is a growing
image, thus increasing the focused region and improving image sharp­ concern that the Zernike fitting method itself may be inaccurate under
ness [5]. This technique allows scanning 25 radial sections of the cornea abnormal conditions. Furthermore, it is very difficult to assess a priori

* Corresponding author at: Department of Structures, Construction and Graphical Expression, Technical University of Cartagena, Campus Muralla del Mar, Car­
tagena 30203, Spain.
E-mail address: [email protected] (F. Cavas).

https://doi.org/10.1016/j.measurement.2024.115613
Received 26 March 2024; Received in revised form 12 August 2024; Accepted 27 August 2024
Available online 30 August 2024
0263-2241/© 2024 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-
nc/4.0/).
A. Ballesta et al. Measurement 240 (2025) 115613

how many terms are necessary to achieve acceptable accuracy in the 0.5 mm [33]. However, these empirical findings do not allow estimating
Zernike reconstruction of any given corneal shape. It is known that to what extent extrapolation/interpolation can be successfully
limiting Zernike analysis to only a few orders may cause incorrect performed.
assessment of the severity of the more advanced stages of keratoconus One possible solution for missing data in computed tomography is to
[9,10]. obtain them from prior information about the surface, which can be
Therefore, various techniques have been developed in recent years to acquired through statistics of the general population [34]. In this way, a
reconstruct and complete elevation and corneal thickness maps from statistical interpretation of the surface restoration problem in computed
incomplete or noisy data. tomography can be reached.
For corneal reconstruction, zonal or modal methods can be utilized. Therefore, in this study, the use of data extrapolation from popula­
In zonal methods, the data collection domain is divided into smaller tion data is proposed to complete elevation/thickness corneal maps and
subdomains, such as triangles. Then, the surface is independently thus improve digital models. It is expected that this will address some of
approximated in each subdomain using splines, particularly B-splines, the limitations found in existing methods, improving their accuracy, and
known for their numerical stability [11,12]. These methods are thus increasing the range of visualization in those areas where the Sirius
commonly used in computer-aided geometric design. This approach has tomograph shows erroneous data (Fig. 1).
been followed by our research group in previous works [13,14]. Other It is expected that the results presented will have significant impli­
authors, such as Chen et al. [15], have applied Deep Learning techniques cations in the field of ophthalmology, helping to improve the accuracy of
along with the zonal approach, presenting a deep learning architecture digital models used for the evaluation and treatment of eye conditions.
called Mesh Reconstruction Network (MR-Net), which enables accurate
reconstruction of three-dimensional meshes in real-time, even with 2. Patients and methods
missing data and sparse annotations for detecting cardiac pathologies.
However, they do not comment on its degree of sensitivity, and although 2.1. Study population
the authors indicate that it could be used for any tissue, their study fo­
cuses on hearts and has not been applied to corneal tissue yet. This cross-sectional study comprised a total number of 86 anterior
Conversely, modal methods involve approximating the surface as a and 86 posterior corneal surfaces, obtained from 43 eyes of 43 adult
linear combination of functions from a set or dictionary defined by a set patients (21 men, mean age 35.10, 22 women, mean age 33.95), 33 of
of parameters. Decisions must be made regarding the selection of them with a confirmed diagnosis of KC (10 Grade I, 10 Grade II, 4 Grace
functions to use, the value of their parameters, and the number of III and 9 Grade IV) and 10 controls with normal corneal topography.
functions needed to recover relevant information without overfitting All subjects underwent ophthalmological examinations including
(the model selection problem). CDVA assessment, slit-lamp bio-microscopy, applanation tonometry
While zonal methods offer flexibility and precision in data fitting, (Goldman), dilated fundus examination and retinoscopy. Participants
they require substantially greater computational resources and encode were requested to remove their contact lenses prior to the measurements
the final shape in a larger amount of data, thus lacking the simplicity of for at least two weeks (for soft contact lenses) or three weeks (for hard
the modal approach, which provides valid functional expressions contact lenses). A single experienced optometrist performed the topog­
throughout the domain, suitable for subsequent calculations such as ray raphy measurements (Sirius System®, CSO, Florence, Italy). Tests with
tracing. Some examples of modal methods include Zernike polynomials the best acquisition quality (green-colored checkmark) were used for the
[16], pseudo-Zernike polynomials [17], Bessel functions [18], Hermite- study.
Gauss polynomials [19], and orthogonal Fourier-Mellin polynomials All patients provided written consent to participate in the study, and
[20]. it was conducted in accordance with the guidelines set forth by the
There are various interpolation or fitting techniques [21–27] that Helsinki Declaration (7th edition, Fortaleza 2013), and was approved by
allow reconstructing a complete surface, such as the cornea, from the ethics committee of the Polytechnic University of Cartagena
incomplete and noisy measurements. Among them, one of the most used (CEI21_001).
is fitting using Zernike polynomials [28,29] due to their simplicity and
their relationship with corneal aberrations. However, there is evidence 2.2. KC group
suggesting that these techniques may not be effective in capturing
important topographic features in irregular corneas [30]. This evidence A single experienced cornea specialist (JLA) verified the diagnosis of
was studied by Baraya et al., who found limitations in the reconstruction KC based on the combination of the following findings; retinoscopic and
of corneal tomography in rabbits using Zernike polynomials with a biomicroscopic signs of KC (such as scissoring, Vogt’s striae, Fleischer’s
Pentacam device, showing that they do not perform well for orders ring, Munson’s sign and Rizzuti’s phenomenon), presence of typical
higher than 12 for the anterior surface and 10 for the posterior surface topographical patterns for KC on axial curvature map (round, oval, su­
[10], which is consistent with the results published by Wei et al. [31]. perior steep, inferior steep, irregular, inferior-steep asymmetric bowtie,
Martinez-Finkelshtein et al. chose to compare different modal recon­ superior-steep asymmetric bowtie and symmetric or asymmetric bowtie
struction algorithms, observing that Zernike polynomials (especially with skewed radial axes (SRAX) > 21 degrees), central/paracentral or
Bhatia-Wolf polynomials) constitute a reliable reconstruction method inferior focal steepening (anterior and/or posterior) and/or corneal
for a surface with non-severe aberrations and a small superficial regu­ thinning, 3-mm inferior-superior (I-S) mean keratometric difference >
larity index (SRI) [9]. However, they fail to capture small deformations 1.4 D. Absence of peripheral KC was also checked. Eyes were graded
of the anterior surface of a synthetic cornea, leading them to propose a regarding disease severity based on the Amsler-Krumeich KC classifi­
new combined approach using an adaptive algorithm, capable of cation system [35].
modeling the cornea from keratometry data capturing smaller details
than with standard procedures [32]. 2.3. Control group
The random and incomplete nature of measurement data can create
an illusion of complete knowledge of the field once the map is restored, The control group included randomly-selected (computer generated)
when in reality important features may be lost or contain measurement one eye of 10 subjects with normal clinical examination and normal
artifacts nonexistent on the actual surface. This uncertainty is especially topography. None of the subjects showed above-mentioned abnormal
relevant in areas where data has not been collected. findings and all had a CDVA≥1.0 Snellen equivalent.
Some analyses of data extrapolation to areas without information Poor compliance to topography measurements, low test quality,
show notable quality in the method if the extrapolation does not exceed previous history of anterior segment surgery, corneal scarring, infection

2
A. Ballesta et al. Measurement 240 (2025) 115613

Fig. 1. Summary of the procedure followed to reach elevation map reconstruction by interpolation methods, using population data.

and any corneal thinning disorders were accepted as the exclusion in which the coefficients represent the free parameters of the model. In
criteria. the case of the anterior surface, since there is no offset (the data starts at
(0,0)), the coefficient C0 is set to 0. The number of points p to be
2.4. Procedure adjusted is generally greater than the number of coefficients of the
polynomial, leading to an overdetermined system (more equations than
To initiate this study on reconstructing missing data in Sirius to­ unknowns), which must be solved so that the resulting polynomial does
pographies from general population statistics, a two-step procedure was not exactly interpolate the data, but rather adjusts them in a sense of
considered (Fig. 2). minimizing the square of the fitting error. It is worth noting that
Firstly, it was deemed necessary to estimate the average of missing although one might consider choosing a high polynomial order to ach­
data per radius across all available eyes in the database (5 8 0), based on ieve a very flexible mathematical model, in practice, it has been shown
each degree of disease, and for both corneal surfaces (anterior and that orders higher than 4 lead to what is known as “overfitting” or
posterior). This aimed to establish the trend of data loss in both the “Runge’s phenomenon” [37], producing unnatural and overly “oscil­
anterior and posterior surfaces, allowing us to determine the optimal lating” surfaces. Therefore, in this case, Grade IV is chosen for the
number of data points to remove based on radius, in order to validate the polynomial fit.
selected reconstruction method for our specific database. Interpolation/extrapolation is carried out from the anterior rings, so
Next, the second step involved setting, for each radius, that exact the fourth-degree polynomial method is applied to perform a partial
number of data points to be removed, denoted as N. corneal reconstruction on these rings. Once the partial corneal recon­
Afterwards, N data points were removed from each ring in a sample struction is done, we will use the coefficients to obtain the interpolation/
of 43 eyes selected per group (anterior cornea group and posterior extrapolation of the heights p(x, y) at each point based on the two-
cornea group) out of the total 580 eyes that make up the general data­ dimensional coordinates.
base. The set of N removed data points was rotated 256 times to create
all possible combinations for each radius. 2.4.1.2. Biharmonic spline interpolation. Biharmonic spline interpola­
It should be noted that all sets of corneal data used were directly tion is a method used to fit a smooth surface to a set of scattered points in
extracted from the tomograhper’s raw datasets, which are unprocessed a three-dimensional space. It is based on constructing an interpolation
data coming from the artificial vision algorithms of Sirius topographer, function that minimizes the biharmonic energy between the given points
and therefore no process was done by any internal Sirius algorithm. [38,39].
Given a scattered set of “n” points with coordinates (xi, yi, zi), the
2.4.1. Interpolation methods function to be evaluated can be expressed as:
Mathematically, the interpolation of point clouds obtained by the ∑n
corneal topographer is a well-known problem that can be approached f(x, y) = i=1
zi ϕi (x, y)
using various methods:
where ϕi (x, y) is a basis function that depends on the location of the
2.4.1.1. Polynomial interpolation. Is one of the most easy and widely control points, and zi the heights associated with them. In this case, a
used fitting techniques [36]. Let p(x, y) be a surface defined by the linear basis function based on triangulation has been chosen, so ϕi (x, y)
function: assigns linear weights to the control points based on the distance be­
tween (x,y) y (xi, yi).
p(x, y) = C0 + C1 • x + C2 • y + C3 • x2 + C4 • x • y + C5 • y2 + C6 Finally, the goal is to minimize the biharmonic energy, which is
• x3 + C7 • x2 • y + C8 • y2 • x + C9 • y3 + C10 • x4 + C11 • x3 defined as the integral of the fourth magnitude of the Laplacian of the
interpolating function:
• y + C12 • x2 • y2 + C13 • x • y3 + C14 • y4

3
A. Ballesta et al. Measurement 240 (2025) 115613

Fig. 2. Graphical description of the procedure.

4
A. Ballesta et al. Measurement 240 (2025) 115613

∫∫
⃒ 2 ⃒
⃒∇ f(x, y) ⃒2 dxdy surfaces (anterior and posterior) was limiting the reconstruction due to
E=
Ω having a higher absence of data (Fig. 3).
It was found that the corneas were complete in all cases up to at least
where Ω is the interpolation domain. This energy minimization can be a radius of 3.6 mm. From this radius onwards, there was a progressive
achieved by various techniques, with the least squares method being increase in data loss, mainly due to erroneous or unmeasured data, with
chosen in this case. Interpolation/extrapolation is carried out from the the frequency and distribution shown in Table 1.
anterior rings, so the biharmonic spline interpolation method is applied Observing the table, it was also determined that starting from a
to perform a partial corneal reconstruction on these rings and obtain the radius of 4.8 mm, erroneous data on the posterior surface are much more
heights of the desired points based on the two-dimensional coordinates. abundant than on the anterior surface, making the posterior surface the
limiting factor when estimating the goodness of the reconstruction.
2.4.1.3. RBF interpolation. Radial basis functions are often used for As a result, based on the results obtained in Table 1, the intervals of
interpolating scattered data obtained through 3D scanning. Methods data to be removed selected were, for the posterior surface:
that use radial basis functions originated to solve the problem of exact [100;120;140;180;200;220;240]; and for the anterior surface:
interpolation in multiple variables. They can generally achieve a smooth [40;50;60;80;100;120;160].
fit of the surface and have a wide variety of basis functions, such as The performance of each of the proposed reconstruction methods
Gaussians, multiquadric inverse, quadratics, inverse multiquadric, or based on the number of eliminated data at each radius was separately
thin-plate splines [40]. studied as a preliminary step of the study, and descriptive representa­
Typically, a parametric surface can be expressed as a sum of indi­ tions can be found in the annex of Supplementary Material Annex 1.
vidual RBFs multiplied by a weighting coefficient [41]: Fig. 4 shows a summary of the results obtained for the anterior and
N
posterior surfaces, for the different degrees of the disease.
∑
S(u, v) = ρi (u, v) • wi The result of applying the selected reconstruction methods is shown
i=1 in Tables 2 and 3 for the anterior and posterior surfaces, respectively.
Fig. 5 shows the result of the reconstruction of the anterior and
The radial basis function technique involves a set of N basis functions, posterior surfaces, as well as the real and reconstructed thickness map,
one for each data point, whose characteristics for the cubic RBF can be for an eye with KC Grade IV.
expressed as: As a mean of evaluating repeatability, a graphical analysis based on
histograms was made. These histograms show the frequency of each
ρi = ‖u − ui ‖3
MAPE% for the fixed set of N eliminated points per radius rotated 256
wi is defined as the weighting coefficients for the RBF. Generally, this times to evaluate it also in positions in which points are known.
coefficient can be n-dimensional and can be considered a control Therefore, each iteration can be considered a new eye, and for each
mechanism since it represents an independent variable that influences radius a total number of N x 256 points are evaluated, making for the
the shape of the final surface. Interpolation/extrapolation is carried out lowest number of points removed (40 in healthy eyes) a total number of
from the anterior rings, so the cubic RBF interpolation method is applied 102.400 points.
to perform a partial corneal reconstruction on these rings and obtain the
heights of the desired points based on the two-dimensional coordinates. 4. Discussion

2.4.2. Surface reconstruction testing Fig. 4 shows how, as the radius increases, and therefore the number
Then, the reconstruction capability of three models (fourth-degree of missing data increases, the error in extrapolation also increases.
polynomial, biharmonic spline, and cubic RBF) was tested on the 43 However, this increase is not linear, at least for all methods, and it is
selected eyes where the previously fixed N data points had been observed that there is no method that outperforms the others for all
removed, evaluating the error incurred in the reconstruction using the degrees of the disease, as those that perform well in diseased eyes give
parameter MAPE (Mean Average Percentage Error), defined as follows worse results for healthy ones. This is particularly evident for the pos­
(eq.1): terior surface.
In Tables 2 and 3, it is observed that, by setting an acceptable MAPE
N ⃒ ⃒
100 ∑ ⃒xi − ̂x i ⃒⃒ % below 5 % for extrapolated radii, the anterior surface can be recon­
MAPE = × ⃒ (1)
N i=1
⃒ x ⃒
i structed from radius 4 to radius 5.2, while the posterior could be
reconstructed from radius 3 to radius 4.2, that is, it is possible to
where, reconstruct up to 6 steps of posterior radii from the starting one, with
sufficient confidence. These results coincide with those shown by some
Xi = actualvalueobtainedfromtheelevationmapofthecornealtopography
authors, who indicate that measurements related to the anterior surface
are more accurate and repeatable than those of the posterior [31], hence
X
̂ i = calculatedvalueobtainedfromtheextrapolationmodels.
a radius further away from the central zone can be reached without
exceeding the 5 % MAPE. The fact that the posterior surface can only be
N=number of calculated points.
reconstructed up to radius 4.2 is consistent with the results observed in
The interpretation of (eq.1) is the average of the percentage errors
other studies and can be explained by multiple reasons, such as optical
produced by the differences between the actual value and the calculated
distortion created by aberrations [20], irregularity of the tear film [31],
value.
or variation in the magnification ratio of the posterior cornea [42],
which cause greater noise in the data [43], and therefore a faster vari­
3. Results
ation of the MAPE% as the distance to the central corneal zone increases.
The three mentioned methods show a very acceptable MAPE% for
The procedure outlined above was put into practice in a statistical
the anterior surface, with MAPE% values below 2, even for the furthest
study based on 580 eyes from the IBERIA keratoconus database (pro­
radii, although a graphical representation of the reconstructed radii
vided by the VISSUM Innovation, Cornea, Cataract and Refractive Sur­
(Fig. 6) allows us to obtain additional conclusions.
gery Unit, Alicante, Spain, which is also affiliated with Miguel
As can be observed, the polynomial method presents continuity is­
Hernandez University). The study aimed to determine up to which
sues, an effect already observed by other authors [44], so it is necessary
radius complete data were available and which of the two corneal
to smooth the contour points, with this method tending to average the

5
A. Ballesta et al. Measurement 240 (2025) 115613

Fig. 3. Scatterplot of the average number of missing data as a function of radius for all G2 eyes present in IBERIA database (N=78) for: a) Anterior corneal surface; b)
posterior corneal surface.

surface. On the other hand, splines effectively represent the shape of the
Table 1
anterior radii, but it is important to note that the accuracy of the
Average of lost data depending on the radius, for both anterior and posterior
biharmonic B-spline approximation is easily influenced by the presence
surfaces.
of noise [45], while the cubic RBF method overpowers the curves,
Average missing Data per radius (Anterior Corneal Surface)
generating less uniform surfaces, but at the same time is less sensitive to
Radius Healthy KC Grade KC Grade KC Grade KC Grade noise [46].
(mm) eyes I II III IV Regarding repeatibility of the method, all histograms show a similar
2.4 0 10 0 0 0 distribution, balanced around the mean value, with all of them showing
2.6 0 14 0 0 5 right skewness (positive skew), meaning there is a long tail towards
2.8 0 15 5 1 10
higher MAPE % values, but with a very low frequency.
3.0 0 22 8 31 26
3.2 0 28 9 36 24 In all cases, most of the data is concentrated in the lower range of
3.4 0 27 19 28 25 MAPE %, indicating that the vast majority of the estimates have rela­
3.6 9 28 21 49 38 tively small errors. Although most of the data is near the mean, there are
3.8 17 27 28 64 33 some values that are farther away, indicating the presence of outliers.
4.0 22 25 28 48 37
4.2 22 29 31 55 41
These outliers most probably correspond to intrinsic measurement er­
4.4 26 30 34 41 45 rors, or typical peripheral measurement error values. A the standard
4.6 29 35 41 50 52 deviation is low and most of the data falls within the mean and two times
4.8 34 40 46 65 58 the standard deviation, it can be said that the results are repeatable. The
5.0 39 45 53 73 69
full set of histograms for all 3 methods has been included as Supple­
5.2 46 51 59 87 78
5.4 65 59 73 94 87 mentary Material Annex 2.
5.6 89 74 95 94 99
5.8 123 98 122 119 123 5. Conclusions
6.0 155 130 149 145 147

In this article, the use of data extrapolation from population-based

Average missing Data per radius (Posterior Corneal Surface) data has been proposed to complete corneal elevation/thickness maps
Radius Healthy KC Grade KC Grade KC Grade KC Grade
and thereby improve digital models. It is expected that this approach
(mm) eyes I II III IV
2.4 0 1 0 0 2
will address some of the limitations found in existing methods,
2.6 0 1 0 0 2 enhancing their accuracy and extending the range of visualization in
2.8 0 1 0 0 3 areas where the Sirius tomograph displays erroneous data.
3.0 0 1 0 3 4 A possible limitation of this study is that the extrapolated version of
3.2 0 2 1 3 5
the data obtained by Sirius was used in those cases in which a part of the
3.4 1 3 2 5 7
3.6 2 4 3 8 10 surface could not be measured in a reliable manner, something that was
3.8 3 6 5 11 14 more frequent for posterior corneal surfaces of patients with KC. This is
4.0 5 9 7 12 20 inevitable, due to the limitations of the measurement technology, and
4.2 10 14 10 14 27
accordingly, a set of cautions need to be taken when interpreting the
4.4 20 22 14 21 41
4.6 39 37 21 29 58
results based on these data. Also, this study only includes elevation data
4.8 67 60 26 38 79 obtained by Sirius tomographer, so it is likely that these results could not
5.0 100 91 34 49 107 be easily applied to data obtained from other topographic systems that
5.2 136 128 42 58 140 are based on other measurement principles and/or use smoothening
5.4 171 164 53 70 171
algorithms.
5.6 197 191 75 86 196
5.8 218 214 103 109 219 Also, it should be noted that this method depends on the database
6.0 237 231 143 145 233 used to fix the number of points to remove in each radius for the vali­
dation method, so although the method would essentially still be valid
for other databases, using a different one could provide different results.
The results presented could have significant implications in the field
of ophthalmology, such as the reconstruction of personalized corneal

6
A. Ballesta et al. Measurement 240 (2025) 115613

Fig. 4. Comparison of interpolation/extrapolation methods. The x-axis shows the radius at which the reconstruction is performed, and the y-axis shows the average
MAPE in % for each N selected per radius.

7
A. Ballesta et al. Measurement 240 (2025) 115613

Table 2
Mean % and SD of the MAPE value for the different reconstruction methods. Anterior surface.
RADIUS (mm)

4.0 4.2 4.4 4.6 4.8 5.0 5.2

ANTERIOR SURFACE %MAPE Healthy Polynomial 0.141 ± 0.237 ± 0.347 ± 0.498 ± 0.668 ± 0.860 ± 0.973 ±
(mean ± SD) 0.139 0.274 0.379 0.528 0.666 0.852 1.064
Splines 0.087 ± 0.214 ± 0.387 ± 0.556 ± 0.803 ± 1.087 ± 1.279 ±
0.064 0.256 0.466 0.731 1.039 1.383 1.737
Cubic 0.192 ± 0.505 ± 0.835 ± 1.099 ± 1.345 ± 1.592 ± 1.639 ±
0.100 0.290 0.526 0.824 1.105 1.419 1.687
KC Grade I Polynomial 0.206 ± 0.318 ± 0.492 ± 0.613 ± 0.735 ± 0.918 ± 1.010 ±
0.208 0.301 0.447 0.614 0.768 0.969 1.185
Splines 0.081 ± 0.216 ± 0.387 ± 0.532 ± 0.728 ± 1.002 ± 1.174 ±
0.068 0.229 0.457 0.697 0.974 1.298 1.621
Cubic 0.185 ± 0.463 ± 0.807 ± 1.063 ± 1.338 ± 1.612 ± 1.676 ±
0.093 0.271 0.518 0.828 1.090 1.383 1.676
KC Grade Polynomial 0.293 ± 0.436 ± 0.593 ± 0.788 ± 0.983 ± 1.194 ± 1.292 ±
II 0.303 0.473 0.633 0.790 0.988 1.242 1.478
Splines 0.108 ± 0.261 ± 0.429 ± 0.648 ± 0.886 ± 1.176 ± 1.368 ±
0.106 0.242 0.398 0.663 0.951 1.275 1.577
Cubic 0.165 ± 0.400 ± 0.693 ± 0.905 ± 1.113 ± 1.315 ± 1.348 ±
0.102 0.275 0.511 0.749 1.003 1.245 1.452
KC Grade Polynomial 0.329 ± 0.466 ± 0.610 ± 0.727 ± 0.842 ± 1.026 ± 1.108 ±
III 0.288 0.419 0.563 0.729 0.922 1.161 1.388
Splines 0.162 ± 0.340 ± 0.592 ± 0.823 ± 1.081 ± 1.422 ± 1.591 ±
0.221 0.363 0.577 0.843 1.139 1.465 1.789
Cubic 0.164 ± 0.345 ± 0.544 ± 0.696 ± 0.870 ± 1.059 ± 1.120 ±
0.173 0.282 0.446 0.636 0.843 1.057 1.262
KC Grade Polynomial 0.529 ± 0.729 ± 0.902 ± 1.076 ± 1.305 ± 1.527 ± 1.601 ±
IV 0.513 0.704 0.899 1.108 1.366 1.649 1.932
Splines 0.200 ± 0.447 ± 0.754 ± 1.049 ± 1.380 ± 1.704 ± 1.882 ±
0.269 0.493 0.776 1.084 1.435 1.823 2.218
Cubic 0.148 ± 0.313 ± 0.485 ± 0.628 ± 0.814 ± 1.028 ± 1.125 ±
0.209 0.344 0.506 0.669 0.874 1.114 1.333

Table 3
Mean % and SD of the MAPE value for the different reconstruction methods. Posterior surface.
RADIUS (mm)

3.0 3.2 3.4 3.6 3.8 4.0 4.2

POSTERIOR SURFACE %MAPE Healthy Polynomial 0.131 ± 0.230 ± 0.330 ± 0.405 ± 0.545 ± 0.724 ± 0.966 ±
(mean ± SD) 0.107 0.188 0.282 0.387 0.519 0.698 0.934
Splines 0.600 ± 1.401 ± 2.230 ± 2.729 ± 3.473 ± 4.142 ± 4.711 ±
0.175 0.585 1.141 1.877 2.475 3.083 3.680
Cubic 0.299 ± 0.777 ± 1.362 ± 1.842 ± 2.571 ± 3.334 ± 4.089 ±
0.135 0.398 0.779 1.303 1.807 2.374 2.996
KC Grade I Polynomial 0.225 ± 0.420 ± 0.637 ± 0.771 ± 0.992 ± 1.230 ± 1.494 ±
0.200 0.419 0.704 0.985 1.298 1.639 2.023
Splines 0.548 ± 1.302 ± 2.093 ± 2.560 ± 3.249 ± 3.862 ± 4.380 ±
0.179 0.564 1.094 1.787 2.345 2.904 3.445
Cubic 0.289 ± 0.780 ± 1.395 ± 1.893 ± 2.643 ± 3.427 ± 4.209 ±
0.139 0.401 0.793 1.329 1.834 2.393 2.995
KC Grade Polynomial 0.396 ± 0.680 ± 0.968 ± 1.120 ± 1.403 ± 1.697 ± 2.020 ±
II 0.386 0.644 0.959 1.271 1.630 2.039 2.513
Splines 0.488 ± 1.116 ± 1.756 ± 2.125 ± 2.663 ± 3.125 ± 3.498 ±
0.230 0.650 1.162 1.746 2.258 2.755 3.217
Cubic 0.294 ± 0.753 ± 1.320 ± 1.785 ± 2.481 ± 3.204 ± 3.919 ±
0.174 0.490 0.892 1.390 1.885 2.423 2.995
KC Grade Polynomial 0.349 ± 0.593 ± 0.821 ± 0.982 ± 1.283 ± 1.608 ± 1.983 ±
III 0.436 0.677 0.926 1.232 1.623 2.039 2.478
Splines 0.395 ± 0.952 ± 1.431 ± 1.698 ± 2.120 ± 2.488 ± 2.793 ±
0.275 0.857 1.202 1.657 2.105 2.531 2.908
Cubic 0.288 ± 0.751 ± 1.203 ± 1.555 ± 2.111 ± 2.662 ± 3.180 ±
0.215 0.626 0.935 1.420 1.910 2.421 2.933
KC Grade Polynomial 0.474 ± 0.841 ± 1.271 ± 1.573 ± 2.097 ± 2.713 ± 3.438 ±
IV 0,417 0.789 1.305 1.833 2.454 3.191 4.071
Splines 0.310 ± 0.704 ± 1.101 ± 1.345 ± 1.718 ± 2.064 ± 2.367 ±
0.208 0.542 0.942 1.352 1.739 2.101 2.429
Cubic 0.173 ± 0.443 ± 0.765 ± 1.043 ± 1.474 ± 1.939 ± 2.412 ±
0.135 0.361 0.658 0.994 1.353 1.739 2.139

models, that can be later used in applications such as personalized or improvements in the diagnosis and follow-up of diseased eyes that
contact lenses fabrication, biomechanical studies that may better were previously not modellable due to the lack of data in the periphery.
simulate incisions in corneal tissue, leading to better surgical practices;

8
A. Ballesta et al. Measurement 240 (2025) 115613

Fig. 5. Result of the reconstruction of the anterior and posterior surfaces, as well as the real and reconstructed thickness maps, for an eye with KC Grade IV.

Fig. 6. Visual comparison of the reconstructed radii by the different interpolation/extrapolation methods.

9
A. Ballesta et al. Measurement 240 (2025) 115613

Funding [13] J.S. Velázquez, F. Cavas, D.P. Piñero, F.J.F. Cañavate, J.A. Del Barrio, J.L. Alio,
Morphogeometric analysis for characterization of keratoconus considering the
spatial localization and projection of apex and minimum corneal thickness point,
This publication was carried out within the framework of the project J. Adv. Res. 24 (2020) 261–271.
“Desarrollo y validación de un nuevo concepto de caracterización bio­ [14] F. Cavas-Martínez, D.G. Fernández-Pacheco, E. De la Cruz-Sánchez, J. Nieto
mecánica-morfofuncional de la córnea” (Reference No. DTS21/00103). Martínez, F.J. Fernández Cañavate, A. Vega-Estrada, et al., Geometrical custom
modeling of human cornea in vivo and its use for the diagnosis of corneal ectasia,
This project was funded by Instituto de Salud Carlos III (ISCIII) and PLoS One 9 (10) (2014) e110249.
cofunded by the European Union. This study was also carried out in [15] X. Chen, N. Ravikumar, Y. Xia, R. Attar, A. Diaz-Pinto, S.K. Piechnik, et al., Shape
collaboration with the Association for the Integration of the Disabled in registration with learned deformations for 3D shape reconstruction from sparse
and incomplete point clouds, Med. Image Anal. 74 (2021) 102228.
the Comarca del Mar Menor (AIDEMAR, collaboration protocol UPCT- [16] H. Lyu, Y. Huang, B. Sheng, Z. Ni, Absolute optical flatness testing by surface shape
AIDEMAR). reconstruction using Zernike polynomials, Opt. Eng. 57 (9) (2018) 094103.
[17] J. Schwiegerling, Scaling pseudo-Zernike expansion coefficients to different pupil
sizes, Opt. Lett. 36 (16) (2011) 3076–3078.
CRediT authorship contribution statement [18] J.C. Carr, R.K. Beatson, B.C. McCallum, W.R. Fright, T.J. McLennan, T.J. Mitchell,
editors. Smooth surface reconstruction from noisy range data2003.
Alejandro Ballesta: Writing – original draft, Validation, Software, [19] A. Wünsche, General Hermite and Laguerre two-dimensional polynomials, J. Phys.
A Math. Gen. 33 (8) (2000) 1603.
Methodology, Formal analysis, Data curation, Conceptualization. José [20] J. Wang, X. Li, Z. Wang, P.G. Davey, Y. Li, L. Yang, et al., Accuracy and reliability
S. Velázquez: Writing – review & editing, Writing – original draft, of orthogonal polynomials in representing corneal topography, Medicine in Novel
Visualization, Validation, Supervision, Data curation, Conceptualiza­ Technology and Devices. 15 (2022) 100133.
[21] Schaap WE, Van De Weygaert R. Continuous fields and discrete samples:
tion. Jorge L. Alió del Barrio: Writing – review & editing, Writing –
reconstruction through Delaunay tessellations. arXiv preprint astro-ph/0011007.
original draft, Resources, Methodology, Investigation, Formal analysis. 2000.
Francisco Cavas: Writing – review & editing, Writing – original draft, [22] L. Wang, D. Chernyak, D. Yeh, D.D. Koch, Fitting behaviors of Fourier transform
and Zernike polynomials, J Cataract Refract Surg 33 (6) (2007) 999–1004.
Supervision, Project administration, Investigation, Funding acquisition,
[23] D.R. Iskander, M.R. Morelande, M.J. Collins, B. Davis, Modeling of corneal surfaces
Conceptualization. with radial polynomials, IEEE Trans. Biomed. Eng. 49 (4) (2002) 320–328.
[24] R. Navarro, L. González, J.L. Hernández, Optics of the average normal cornea from
general and canonical representations of its surface topography, JOSA a. 23 (2)
Declaration of competing interest
(2006) 219–232.
[25] T. Gu, H. Lin, D. Tang, S. Lin, T. Luo, Curve and surface reconstruction based on
The authors declare that they have no known competing financial MTLS algorithm combined with k-means clustering, Measurement 182 (2021)
109737.
interests or personal relationships that could have appeared to influence
[26] S. Zhang, Q.X. Chen, J.W. Luo, J. Bai, A fast surface reconstruction method for
the work reported in this paper. fluorescence molecular tomography based on cross-beam edge back projection,
Measurement 46 (4) (2013) 1565–1571.
Data availability [27] Zhang Y, Kim NH, Park C, Haftka RT, editors. Function extrapolation of noisy data
using converging lines2016.
[28] S.D. Klyce, M.D. Karon, M.K. Smolek, Advantages and disadvantages of the Zernike
Data will be made available on request. expansion for representing wave aberration of the normal and aberrated eye,
J. Refract. Surg. 20 (5) (2004).
[29] J. Schwiegerling, Modal reconstruction methods with Zernike polynomials. Slack
Appendix A. Supplementary material Incorporated Thorofare, NJ; 2005. p. S552-S7.
[30] M.K. Smolek, S.D. Klyce, Zernike polynomial fitting fails to represent all visually
Supplementary data to this article can be found online at https://doi. significant corneal aberrations, Invest. Ophthalmol. Vis. Sci. 44 (11) (2003)
4676–4681.
org/10.1016/j.measurement.2024.115613. [31] Y. Wei, B.T. Lopes, A. Eliasy, R. Wu, A. Fathy, A. Elsheikh, et al., Performance of
Zernike polynomials in reconstructing raw-elevation data captured by Pentacam
References HR, Medmont E300 and Eye Surface Profiler, Heliyon. 7 (12) (2021) e08623.
[32] A. Martınez-Finkelshtein, D. Ramos-López, G.M. Castro-Luna, J.L. Alió. An
adaptive algorithm for the cornea modeling from keratometric data. arXiv preprint
[1] Y. Li, D.M. Meisler, M. Tang, A.T.H. Lu, V. Thakrar, B.J. Reiser, et al., Keratoconus
arXiv:10091244. 2010.
diagnosis with optical coherence tomography pachymetry mapping,
[33] G.-m. Dai, Validity of scaling zernike coefficients to a larger diameter for refractive
Ophthalmology 115 (12) (2008) 2159–2166.
surgery, J. Refract. Surg. 27 (11) (2011) 837–841.
[2] R. Ambrósio Jr, A.L.C. Caiado, F.P. Guerra, R. Louzada, A.S. Roy, A. Luz, et al.,
[34] A. Fabrikant, Objective Estimation for Uncertainty of Restoring Corneal
Novel pachymetric parameters based on corneal tomography for diagnosing
Topography Surface, Eye Contact Lens 41 (2) (2015) 121–126.
keratoconus, J. Refract. Surg. 27 (10) (2011) 753–758.
[35] M. Amsler, The“ forme fruste” of keratoconus, Wien. Klin. Wochenschr. 73 (1961)
[3] S. Binder, C.I. Falkner-Radler, C. Hauger, H. Matz, C.G. Glittenberg, Clinical
842–843.
applications of intrasurgical SD-optical coherence tomography, Invest.
[36] Surface fitting | Scilab.
Ophthalmol. Vis. Sci. 51 (13) (2010) 268.
[37] B. Fornberg, J. Zuev, The Runge phenomenon and spatially variable shape
[4] F. Bao, G. Savini, B. Shu, S. Zhu, R. Gao, G. Dang, et al., Repeatability,
parameters in RBF interpolation, Comput. Math. Appl. 54 (3) (2007) 379–398.
reproducibility, and agreement of two Scheimpflug-Placido anterior corneal
[38] D.T. Sandwell, Biharmonic spline interpolation of GEOS-3 and SEASAT altimeter
analyzers for posterior corneal surface measurement, J. Refract. Surg. 33 (8)
data, Geophys. Res. Lett. 14 (2) (1987) 139–142.
(2017) 524–530.
[39] F. Hou, Q. Gao, Y. Song, Z. Wang, Z. Bai, Y. Yang, et al., Deep feature pyramid
[5] F. Cavas-Martínez, S.E. De la Cruz, J. Nieto Martínez, F.J. Fernández Cañavate, D.
network for EEG emotion recognition, Measurement 201 (2022) 111724.
G. Fernández-Pacheco, Corneal topography in keratoconus: state of the art, Eye and
[40] H. Rocha, On the selection of the most adequate radial basis function, App. Math.
Vision. 3 (2016) 1–12.
Model. 33 (3) (2009) 1573–1583.
[6] H.M. Merklinger, Focusing the view camera, Seaboard Printing Limited. 5 (1996).
[41] I. Marinić-Kragić, S. Perišić, D. Vučina, M. Ćurković, Superimposed RBF and B-
[7] Pinero DP, editor Technologies for anatomical and geometric characterization of
spline parametric surface for reverse engineering applications, Integr. Comput.-
the corneal structure and anterior segment: a review2015: Taylor & Francis.
Aided Eng. 27 (1) (2020) 17–35.
[8] M.R. Hee, Artifacts in optical coherence tomography topographic maps, Am. J.
[42] X. Zheng, W. Yang, L. Huang, J. Wang, S. Cao, B. Geraghty, et al., Evaluating the
Ophthalmol. 139 (1) (2005) 154–155.
repeatability of corneal elevation through calculating the misalignment between
[9] A. Martinez-Finkelshtein, A.M. Delgado, G.M. Castro, A. Zarzo, J.L. Alió,
Successive topography measurements during the follow up of LASIK, Sci. Rep. 7 (1)
Comparative analysis of some modal reconstruction methods of the shape of the
(2017) 3122.
cornea from corneal elevation data, Invest. Ophthalmol. Vis. Sci. 50 (12) (2009)
[43] S. Schröder, S. Mäurer, T. Eppig, B. Seitz, K. Rubly, A. Langenbucher, Comparison
5639–5645.
of corneal tomography: repeatability, precision, misalignment, mean elevation,
[10] M. Baraya, J. Moore, B.T. Lopes, R. Wu, F. Bao, X. Zheng, et al., Limitations of
and mean pachymetry, Curr. Eye Res. 43 (6) (2018) 709–716.
reconstructing pentacam rabbit corneal tomography by zernike polynomials,
Bioengineering 10 (1) (2022) 39.
[11] M.A. Halstead, B.A. Barsky, S.A. Klein, R.B. Mandell, A spline surface algorithm for
reconstruction of corneal topography from a videokeratographic reflection pattern,
Optom. Vis. Sci. 72 (11) (1995) 821–827.
[12] B. Büchmann, Accuracy and stability of a set of free-surface time-domain boundary
element models based on B-splines, Int. J. Numer. Meth. Fluids 33 (1) (2000)
125–155.

10
A. Ballesta et al. Measurement 240 (2025) 115613

[44] C. Mineo, S.G. Pierce, R. Summan, Novel algorithms for 3D surface point cloud [46] L.K. Jie, A. Ramli, A.A. Majid, A comparison study between b-spline surface fitting
boundary detection and edge reconstruction, J. Comput. Des. Eng. 6 (1) (2019) and radial basis function surface fitting on scattered points, Jurnal Teknologi. 78
81–91. (6–5) (2016) 103–107.
[45] K.J. Liew, A. Ramli, M.A. Abd, B-spline surface fitting on scattered points, Appl.
Math. Inf. Sci. 10 (1) (2016) 273.

11