Study and characterization of morphogeometric parameters to assist diagnosis of keratoconus

Background In case of significant imperfections on the cornea, data acquisition is difficult and a significant level of missing data could require the interpolation of important areas of the cornea, resulting in a very ambiguous model. The development of methods to define in vivo customised geometric properties of the cornea based only on real raw data is extremely useful to diagnose and assess the progression of diseases directly related to the corneal architecture. The present work tries to improve the prognostic of corneal ectasia creating a 3D customised model of the cornea and analysing different geometric variables from this model to determine which variables or combination of them could be defined as an indicator of susceptibility to develop keratoconus. Methods A corneal geometric reconstruction was performed using zonal functions and retrospective Scheimpflug tomography data from 187 eyes of 187 patients. Morphology of healthy and keratoconic corneas was characterized by means of geometric variables. The performance of these variables as predictors of a new geometric marker was assessed and their correlations were analysed. Results The more representative variable to classify the corneal anomalies related to keratoconus was posterior apex deviation (area under receiver operating characteristic curve > 0.899; p < 0.0001). However, the strongest correlations in both healthy and pathological corneas were provided by the metrics directly related to the thickness, as deviations of the anterior/posterior minimum thickness points. Conclusions The presented morphogeometric approach based on the analysis and custom geometric modelling of the cornea demonstrates to be useful for the characterization and diagnosis of keratoconus disease, stating that geometrical deformation is an effective marker of the ectatic disease’s progression.


Background
Keratoconus is an ectatic corneal disorder characterized by impaired vision and a poor quality of life. It is one of the leading indications for corneal transplantation. This pathology is characterised by a progressive corneal thinning and a structural weakening, resulting in corneal protrusion, irregular astigmatism and a gradual deterioration of the visual performance related to morphology changes in the corneal architecture. Several diagnostic criteria have been defined using a huge variety of techniques and technologies such as the classical keratoconus biomicroscopic signs, the conical protrusion and the infero-superior asymmetry [1].
As it has been demonstrated that small variations of the corneal morphology could induce important changes in the quality of life of the patients [2], recently other parameters such as geometric parameters of the keratoconic corneal surfaces have been analysed [3].
The geometric reconstruction of the corneal surface has experienced significant progress in recent years with the development of new technologies. In clinical practice the ophthalmologists use corneal topographers based on the Scheimpflug technology [4]. This diagnostic equipment gives a matrix of discrete points for the anterior and posterior surface of the cornea [5]. These raw data are not interpolated and are used to generate a geometric model with the modal methods called Zernike polynomials. These Zernike polynomials are defined for all discrete points of the anterior and posterior corneal surfaces for their reconstruction. However, in cases of pathological corneas, as for instance very aberrated corneas, it has been demonstrated that it is very difficult to define the Zernike polynomial order required to get the most relevant information about the corneal surfaces [5][6][7].
An alternative to these geometric models are the zonal methods called B-Spline functions [8]. These B-Spline functions divide the general raw data area used for the surface reconstruction into small more elemental subareas to get more flexible and more precise adjustment of the geometric model. These methods are widely used in Computer Aided Geometric Design (CAGD) for the virtual reconstruction of complex geometrical surfaces and its posterior analysis in different industrial fields [9,10]. In Bioengineering, the zonal methods are used to generate virtual 3D models of the biological structures for different applications [11][12][13][14][15][16]. However, this technology is not extended to the Ophthalmology field, and it has only been used to generate customised models of the corneal biomechanics [3,[17][18][19][20]. In these cases, the raw data obtained were incomplete in the periphery and the limbus areas due to extrinsic errors during the data acquisition process, so the authors had to interpolate data in order to obtain a complete geometric model.
In the case of significant imperfections on the cornea, the data acquisition is difficult and a significant level of missing data could require the interpolation of important areas of the cornea resulting in a very ambiguous model. The development of methods to define in vivo customised geometric properties of the cornea based only on real raw data and without any kind of interpolation is extremely useful in clinical practice for the modelisation of pathological corneas [5].
The present work tries to improve the prognostic of corneal ectasia using a new concept of diagnosis based on a geometric modelisation of the cornea. The Graphical Bioengineering technique used in this study will create a 3D customised model of the cornea and analyse different geometric variables of this model to determine, in a largescale computational trial, which geometric variables or combination of variables could be defined as an indicator of susceptibility to develop keratoconus (Fig. 1).
A complete geometric analysis of the cornea of a huge number of patients has been carried out on raw data using Computer Aided Geometric Design with zonal methods and without any kind of interpolation. Different geometric variables have been defined for the 3D model of the cornea of each patient. These geometric variables have been analysed and compared according to the evolution of the pathology to define the most adequate one for the detection and prognostic of the evolution of the pathology.
The definition and analysis of the geometric profile of each cornea derived from the corneal pathology will enable to define the first changes related to the development of keratoconus and the evolution of the geometric profile according to the keratoconus severity, the corresponding loss of visual quality and the therapeutic treatments indicated.

Methods
This was a retrospective study including 187 eyes of 187 patients ranging in age between 7 and 73 years old. Only one eye from each patient was randomly selected to avoid the interference in the analysis of the correlation that often exists between the two eyes of the same person. This study was conducted at Vissum Corporation in Alicante (Spain). Two groups of eyes were differentiated depending if the keratoconus disease was present or not: control group, including 124 healthy eyes, and keratoconus group, including 63 eyes with the diagnosis of keratoconus. The inclusion criterion for the control group was healthy eyes that did not meet the exclusion criteria and diagnosis according to the standard criteria for keratoconus diagnosis in the keratoconus group [21,22], which is the presence of an asymmetric bowtie pattern in corneal topography, a value of 100 or higher of the KISA index, a central keratometry (K-value) with different cut-off values to keratoconus suspect (> 47.2 D), an inferior-superior asymmetry (I-S value) with a cutoff value of 1.4 D (difference between average inferior and superior corneal powers at 3 mm from the centre of the cornea), as well as other topographic indexes (SRAX, KSS, KPI, CLMI) and at least one keratoconus sign on slit-lamp examination, such as stromal thining, conical protusion on the cornea at the apex, Fleischer ring, Vogt striae or anterior stromal scar. Exclusion criteria in both groups were previous ocular surgery and any other active ocular disease. This study was approved by the Vissum Corporation ethics committee and was performed in accordance with the ethical standards laid down in the Declaration of Helsinki (Seventh revision, October 2013, Fortaleza, Brasil).

Examination protocol
All patients underwent a complete eye examination including the following tests ( Fig. 1): measurement of uncorrected (UDVA) and corrected distance visual acuity (CDVA), manifest refraction, slit-lamp biomicroscopy, PIO, biometrics (axial length, spherical equivalent, minimal corneal thickness, central corneal thickness, anterior chamber depth, white-to-white corneal diameter) and corneal topographic analysis. During this protocol, the Sirius system ® (CSO, Florence, Italy) was used, which is a noninvasive system for measuring and characterizing the anterior segment using a rotating Scheimpflug camera that generates images in three dimensions, with a dot matrix fine-meshed in the center due to the rotation. The images taken during the examination are digitalized in the main unit and transferred to a computer to be analyzed in detail. Gathered Sirius corneal topographies (data from other topographers such as Pentacam can also be handled [23]) are represented as discrete and finite set of spatial points (point cloud surfaces) in the form of two 31 × 256 matrices. Both matrices contain the polar coordinates representative of the anterior and posterior corneal surfaces. These data, used only in the first stage of the topographic data acquisition procedure, are called raw data [5,23]. This warrantied that data were not interpolated or manipulated [23,24], avoiding any proprietary reconstruction software from topographer's manufacturer. All corneal topography files were exported in CSV format. Likewise, all cases were classified according to the Amsler-Krumeich grading system. All measurements were performed by the same experienced optometrists, performing three consecutive measurements and taking average values for posterior analysis.

Diagnosis procedure
The procedure based on zonal methods and proposed in this article consists of two main stages: i) a geometrical modelling stage where the raw data provided by the corneal topographer is used to reconstruct a 3D geometrical model of the cornea using computational geometry techniques, and ii) a geometrical analysis stage where determined geometric variables are extracted from the model and analyzed to characterize the cornea.
1. Extraction of the point clouds from the corneal topographer. The Sirius device used for this study provides two 3D point clouds that conform both the anterior and posterior corneal surfaces, respectively. However, due to Sirius topographer only provides spatial points data in Cartesian format for the anterior corneal surface (not for the posterior surface), data in polar format had to be exported to obtain data from both corneal surfaces. This data was given as a CSV (Comma Separated Values) table where every row represented a circle in the corneal map and every column represented a semi-meridian, giving 256 points for each radius. This way, each i-th row sampled a circle of i*0.2 mm radius on the map, and each j-th column sampled a semi-meridian in the direction of j*360/256° on the map, so each Z value of the matrix [i, j] represented the point P (i*0.2, j*360/256°) in polar coordinates. In order to perform these calculations, exported data were further formatted in Cartesian coordinates by the aid of an algorithm programmed in Matlab ® . 2. Generation and positioning of both corneal surfaces. The two point clouds representing the corneal geometry were imported into the surface reconstruction software Rhinoceros ® v5.0, which uses a mathematical model to generate surfaces based on non-uniform rational B-spline (NURBS) with high accuracy. The surfaces that best fit the point clouds were generated with the Rhinoceros' patch surface function, a reconstruction software option that fits a surface through given curves, meshes, point objects, and point clouds. For this research, this function tried to minimize the nominal distance between the 3D point clouds and the solution surfaces. For this objective, the function was configured by setting the sample point spacing at 256 (number of points for each data ring), the surface span planes at 255 for both u and v directions (the maximum number of span planes that the software permitted), and the stiffness of the solution surface at 10 −3 (mm). This last parameter provides information on how much the best fit plane can be deformed in order to match the input points. This deviation can be calculated later by the software, providing a mean value of the distance error for the solution surface. During this study, an average distance error between the solution surface and the 3D point cloud of about 3.60 × 10 −4 ± 6.43 × 10 −4 mm was obtained. Using this procedure, anterior and posterior corneal surfaces were generated and engaged by their geometrical center and Z axis. 3. Solid modeling. After generation and positioning of both corneal surfaces, a third peripheral surface (the bonding surface between both sides in the Z-axis direction) was generated and then joined together to form a single surface. The surface reconstructed was then exported to the solid modeling software SolidWorks V 2013 (Dassault Systèmes, Vélizy-Villacoublay, France) to generate a solid model that is representative of the custom and actual geometry of each cornea.
A full process for the geometric reconstruction of a patient-specific cornea that comprises the three mentioned stages takes around 2-3 min.

Second stage: geometrical analysis
The resulting 3D solid model of the cornea was then used to perform an analysis of determined geometric variables (Fig. 3) that are representative of the corneal morphology (Table 1). These variables were later statistically analyzed in order to characterize the cornea.

Statistical analysis
A Kolmogorov-Smirnov test was run to assess the data engagement scores. According to this test and thereafter, a Student's t test or U-Mann-Whitney Wilcoxon test was performed (depending on normality), in order to describe differences between normal and keratoconus groups in all the measurements proposed. Additionally, Kruskal-Wallis (K-W) and Effect Size (ES) tests were used to compare differences and to quantify the For all statistical tests, the same level of significance was used (p < 0.05). Correlation coefficients (Pearson or Spearman depending if normality condition could be assumed) were used to assess the correlation between all different parameters. A linear regression was performed to quantify the strength of the correlation (R 2 ) for both groups studied. A ROC curve analysis was performed in order to obtain the accuracy of the measurements. This ROC curve is a graphical plot that illustrates the performance of a binary classifier system as its discrimination threshold is varied. The curve is created by plotting the true positive rate against the false positive rate at various threshold settings. The accuracy of the test depends on how well the test separates the group being tested into those with and without the disease in question. The area under the ROC curve measures the accuracy: an area of 1 represents a perfect test, and an area of 0.5 represents a worthless test. A rough guide for classifying the accuracy of a diagnostic test is the traditional

Results
From a total of 187 patients, this study included 124 healthy eyes that did not present any ocular pathology [25] and constituted by 69 females (55.6%) and 55 males (44.4%) ranged from 7 to 73 years old, and 63 eyes diagnosed with keratoconus in several grades [26] (53.9% in stage I, 31.7% in stage II, and 14.4% in the most extreme stages, III and IV) and formed by 34 females (53.9%) and 29 males (46.1%) ranged from 14 to 69 years old. Table 2 shows the visual, refractive and morphological outcomes in the two groups of eyes analyzed during the study. No statistically significant differences between groups were found in axial length (p = 0.33, Student's t-test), anterior chamber depth (p = 0.29, Mann-Whitney test), white to white corneal diameter (p = 0.71, Mann-Whitney test) and age (p = 0.09, Mann-Whitney test). Significant differences between groups were On the other hand, all of the morphogeometric variables showed differences between normal and keratoconic eyes, as seen in Table 3. Total corneal volume presents higher values in healthy eyes (p < 0.0001), while anterior and posterior corneal surface areas are lower in the same subjects (p < 0.0001). This pattern of difference can be seen for most of the variables studied: healthy corneas have anterior and posterior apex deviations lower  than keratoconic corneas (p < 0.0001), as occurs with the anterior and posterior minimum thickness point deviations.
Outcomes according to keratoconus severity are shown in Table 4 (clinical features) and Table 5 (morphogeometric parameters), where comparisons are established according to the AK grading system. Additionally, note that calculated effect sizes for each disease stage allows quantifying the degree of change, which is higher for stages III and IV in all of the variables, becoming more evident with the progress of the disease. Table 6 summarizes the statistically significant correlations between all the modeled morphogeometric variables for the normal group, and Table 7 shows the significant correlations for the KC group. In the healthy eyes group two strong correlations (above 0.9) have been achieved: on one hand, between total corneal volume and total corneal surface area (r = 0.953; p = 0.000), and on the other hand between anterior and posterior minimum thickness point deviations (r = 0.996; p = 0.000).
In the same way, two strong correlations (above 0.9) have also been found in the keratoconic eyes group for several analysed morphogeometric variables: between anterior and posterior corneal surface areas (r = 0.921; p = 0.000), and between anterior and posterior minimum thickness point deviations (r = 0.999; p = 0.000).
The predictive value of the modeled variables was established by a ROC analysis (Fig. 4). From the several geometric variables analyzed during the study, the variables that achieved the best results in the diagnosis of the disease with an area under the ROC curve (AUROC) above 0.7 were the following four: anterior corneal surface area (Fig. 3a) (area: 0.853, p < 0.0001, std. error: 0.040, 95% CI 0.762-0.919), with a cutoff value of 43.07 mm 2 , and an associated sensitivity and specificity of 90.27% and 60.01%, respectively; the posterior corneal surface area (Fig. 3a) (area: 0.813, p < 0.0001, std. error: 0.039, 95% CI 0.719-0.891), with a cutoff value of 44.18 mm 2 , and an associated sensitivity and specificity of 91.08% and 44.17%, respectively; anterior apex deviation (Fig. 3b) (area: 0.742, p < 0.0001, std. error: 0.059, 95% CI 0.641-0.875), with a cutoff value of Thus, according to the area under the curve variable calculated for the analysed parameters, it was concluded that the parameter that provides a higher rate of discrimination between normal corneas and corneas with keratoconus is Posterior apex deviation. Nevertheless, there are other relevant statistical differences between healthy and diseased eyes, and most of variables studied differ between groups, making it possible to differentiate with high sensitivity and specificity healthy corneas from those patients diagnosed with keratoconus.

Discussion
The introduction of advanced imaging technologies in clinical practice, such as scanning-slit topography or Scheimpflug tomography has allowed the clinician to perform a more consistent and precise diagnosis of forms of keratoconus, helping to characterize the global geometry of the corneal surface [4]. Several parameters, such as pachymetry, corneal volume or corneal wavefront aberrations, are used in the batteries tests for corneal ectasia diagnosis, providing a characterization of the underlying morphogeometric alteration [21]. However, the geometrical characterization indices proposed by these devices are not easily compatible between different tomographers, generating confusion in the Ophthalmic Community [27][28][29].
On the other hand, up to date, the ability of a specific geometric model to capture diseases on human corneas has not been assessed, either the correlation values between morphogeometric parameters on both normal and pathological groups. Specifically, this computational study provides insight into the complex clinical problem of diagnosing corneal ectatic diseases.
It is well known that the mechanical response of any deformable system is affected by its geometry and material properties. When geometry is fully characterized, it is possible to set up a geometric model of the system, which may be used to analyze the geometric response under original conditions. In this case, conditions are defined by the rupture of the geometric balance due to the existence of a biomechanical weakening, as happens in the keratoconus disease. Keratoconus is a disorder characterized by a progressive thinning of the cornea, which is physically presented in its structure as a protrusion or cone type focal curving that entails a redistribution of its pachymetry and some changes in the anatomical morphology of its surfaces [21]. Other aspect to be taken into account is the geometric characterization based on raw data, which has been previously used by some authors in the Corneal Biomechanics field [20,23] and for diagnosis of corneal diseases [30]. However, these studies resort to data interpolation to obtain a specific model for each patient. Geometric modeling based on raw data that has not been treated by any internal algorithm of the topographer enables an accurate clinical characterization of the human cornea basing on perfectly defined morphological variables in the field of graphic bioengineering.
During this study, two different types of parameters were analysed: (i) clinical features and (ii) morphogeometric variables. During the first analysis, no significant differences in age were found between the different grades of keratoconus (p = 0.31). However, statistically significant differences between the healthy group and the different grades of keratoconus were found for refractive and biometric results (p = 0.0001).
In the second analysis, significant differences were found for the morphogeometric variables between the health and disease groups, as well as between the different evolution grades of the disease. This fact reveals that, even though the curvature radii are smaller on both surfaces in the KC group, the morphogeometric variables register the tendency of the cornea to develop and to maintain its structure in form of meniscus until the most advanced grades of the disease, in which the relationship between both corneal surfaces is significantly modified with a greater increase of the curvature (distance between the corneal apex and the geometric centre of the cornea) of the posterior corneal surface compared to the anterior corneal surface. This trend is on line with the tendency reported by other studies [31], where the posterior-anterior radius of curvature is analysed for groups with healthy eyes and keratoconic eyes according to the Amsler-Krumeich classification. Specifically for keratoconus, in this study a significant increase of the curvature of the posterior corneal surface was observed with respect to the anterior surface (p < 0.01) in the most advanced disease stage. According to these findings, the severity grade of the keratoconus disease seems to be related to the geometry of the anterior and posterior corneal surfaces. This fact could be associated with the alteration of the biomechanical properties existing in a keratoconic cornea, and more concretely in the most advanced grades of the disease. Theoretically, this biomechanical weakness can make the cornea to be subject of deformation due to the intraocular pressure, affecting largely to the posterior surface curvature.
Regarding the ROC analysis, the posterior apex deviation showed the highest discriminant coefficient (AUROC, 0.899). In accordance with a previous publication [32,33], the posterior curvature might influence on the visual function. We hypothesize that the magnitude of the posterior apex deviation represents the best performance to recognize the geometric profile of the KC stage where the visual acuity starts to impair. Furthermore, we also found a satisfactory level of discriminative ability for anterior corneal surface area (AUROC, 0.853), posterior corneal surface area (AUROC, 0.813) and the anterior apex deviation (AUC, 0.742).
This study also assessed the relationship of the geometrical pattern of both anterior and posterior corneal surfaces. Other authors, for a different purpose, found strong statistical correlations between the anterior and posterior shape factors for keratoconic corneas. Regarding the minimum thickness point deviations from both corneal surfaces, significant differences between groups were found. Interestingly, the strongest correlation value yielded in this investigation for keratoconic eyes was between the anterior and posterior minimun thickness point deviations (R 2 = 0.999; p = 0.000).

Conclusions
This study relies on the use of a reduced number of geometrical parameters obtained from modeling tests of the cornea: anterior corneal surface area, posterior corneal surface area, anterior apex deviation and posterior apex deviation. These variables are sufficient to prove that the variability of the geometric response of human corneas is definitely related to diagnosis of the disease. This method is simplified and more integrative than current diagnostic systems, which analyse separately the anterior and posterior surfaces of the cornea. This observation has a quite relevant implication in view of the prediction of the response to refractive surgery, i.e. the knowledge of the sole geometry is enough to feed keratoconus diagnosis.
The main suggestion derived from this study is to give high priority to the development of non-invasive testing methods that are able to provide through inverse analysis the patient-specific parameters of a sufficiently realistic geometric model of the corneal morphology, which can be obtained with the aid of Computer Aided Geometric Design tools.
This method will allow improving the detection and effects of therapeutic methods used for keratoconus and other corneal ectatic diseases such as post-lasik ectasia. Early studies, currently under publication, have demonstrated the effectiveness of this approach in the early detection of subclinical keratoconus. In a close future, thanks to the analyses of the objective data related to the geometric effect of the intracorneal rings implanted, customized nomograms for the implantation of Intra-Corneal Rings will be developed. Later, the analysis of the correlation between the geometric, visual, biomechanical and clinical effects of intracorneal implants in ectatic corneas will allow the development of new therapies and new concepts of corneal implants. The geometric modeling developed will also allow assessing more accurately the outcomes of the corneal crosslinking techniques and its effectiveness in slowing the development of keratoconus.