Corneal Properties of Keratoconus Based on Scheimpflug Light Intensity Distribution

Abstract

Purpose: We introduce a new approach to assess the properties of corneal microstructure in vivo of healthy control and keratoconus eyes, based on statistical modeling of light intensity distribution from Scheimpflug images.

Methods: Twenty participants (10 mild keratoconus and 10 control eyes) were included in this study. Corneal biomechanics was assessed with a commercial Scheimpflug camera technology. Sets of 140 images acquired per measurement were exported for further analysis. For each image, after corneal segmentation, the stromal pixel intensities were statistically modeled, leading to parametric time-series that characterizes distributional changes during the measurement. From those time series, a set of 10 newly introduced parameters (microscopic parameters) was derived to discriminate normal from keratoconic corneas and further compared against clinical parameters available from the same measuring device, including central corneal thickness, IOP, and deformation amplitude (macroscopic parameters).

Results: Biomechanical microscopic parameters extracted from statistical modeling of light intensity distribution were good discriminators between mild keratoconus and control eyes (Mann-Whitney U test, P < 0.05/N [Bonferroni]). The combination of available macroscopic and novel microscopic parameters was the most successful tool to differentiate between keratoconus and control eyes with no misclassifications.

Conclusions: For the first time to our knowledge, a set of parameters related to corneal microstructure, acquired from statistical modeling of light intensity distribution of dynamic Scheimpflug image acquisition was introduced. This novel approach showed the potential of combining macroscopic and microscopic corneal properties derived from a single clinical device to discriminate successfully between mild keratoconus and control eyes.

The cornea is a complex anisotropic composite structure with nonlinear elastic and viscoelastic properties.¹ Viscoelastic corneal properties depend on the internal corneal structure and composition and greatly fluctuate between subjects. Age,^2,3 pregnancy,^4–6 presence of certain hormones,^4,7,8 and illnesses related to collagen synthesis^8,9 are only a few factors that have influenced corneal biomechanics. While the cornea is more commonly thought of as an optical rather than a mechanical system, corneal biomechanics has a critical role in a number of ophthalmic conditions and processes,¹⁰ including keratoconus.^11,12 

Corneal mechanical stability is compromised in keratoconus,¹³ potentially leading to progressive macroscopic morphologic changes. For the sake of better understanding the mechanism underlying keratoconus development and to enhance the rates of early diagnosis, previous works investigated the macroscopic biomechanical differences between healthy and keratoconus eyes.^14–18 In vivo macroscopic measurements of corneal biomechanics are minimally invasive, and can be performed in clinical practice using commercially available air-puff systems,¹¹ such as the Corvis ST (Oculus, Wetzlar, Germany). This noncontact tonometer emits an air-puff and uses a high-speed Scheimpflug camera to record a series of images depicting a cross-section of the cornea during deformation. The deformation profiles then are analyzed to derive several geometry-based parameters.¹⁹ 

A detailed mechanical description of corneal biomechanics should involve an analysis of microscopic and macroscopic parameters, but an analysis of microscopic parameters is not yet clinically available. Recently, it was demonstrated that information about microscopic corneal biomechanical properties can be extracted in vivo through the statistical modeling of optical coherence tomography (OCT) speckle.^20–23 Besides the speckle inherent to OCT measurements, traditionally considered as noise, it was demonstrated that there also is a component that carries information about the microstructure of the tissue.²⁴ Speckle arises from interferometry, the optical technique on which OCT is based. Scheimpflug devices, on the other hand, are not based on interferometry and consequently do not produce speckle. Instead, the light intensity distribution in the corneal Scheimpflug images may have potential in the assessment of microstructural corneal properties in a similar way as OCT.²⁰ 

We introduced a set of parameters related to corneal microstructure, based on the statistical modeling of light intensity distribution of dynamic Scheimpflug images and showed the potential of combining macroscopic and microscopic properties of the cornea derived from a single clinical device to discriminate between keratoconus and control eyes. 

Methods

Subjects and Data Collection

We studied 20 adults (20 eyes; 11 females, nine males) between 18 and 29 years old (mean age ± SD, 23.6 ± 2.9 years). The keratoconus group (KC) was formed by 10 participants, while the other 10 had healthy control eyes with IOP within a narrow range of 15 to 17 mm Hg and corneal astigmatism less than 0.75 diopters (D). Keratoconus patients were recruited from the keratoconus clinic at the Antwerp University Hospital (UZA, Edegem, Belgium). The keratoconic eyes were not crosslinked and were not planned to be crosslinked in the next year at the moment of recruitment, meaning that they represented a rather mild keratoconus. Besides crosslinking, other exclusion criteria included corneal scarring, contact lens wear, known retinal or corneal pathologies (apart from keratoconus), known ocular procedures or treatments, and known systemic diseases (e.g., diabetes, human immunodeficiency virus/acquired immunodeficiency syndrome [HIV/AIDS], or hypertension). In patients who were diagnosed with keratoconus in only one eye, the pathologic eye was selected for measurement. For participants with keratoconus in both eyes and for control subjects, one eye was randomly selected for measurement and further statistical analysis. All participants underwent a comprehensive ophthalmologic examination, including corneal biomechanics with Corvis ST (software version 1.4r1755) and a standardized interview. All assessments were performed at the same time of day during a single visit to decrease the effect of diurnal fluctuation. The study was approved by the Antwerp University Hospital ethical committee and adhered to the tenets of the Declaration of Helsinki. All participants gave written informed consent to participate after the nature and possible consequences of the study were explained. 

Data Analysis

Corvis ST corneal deformation images were obtained every 230 μs and had a fixed size of 200 × 576 pixels, corresponding with an approximate resolution of 15 μm/pixel. After data acquisition, the sets of 140 images acquired per measurement were exported for further analysis. The method of data analysis consisted of two steps: corneal segmentation and statistical modeling of the pixel intensity distribution. Traditional image processing techniques, such as median filter and Canny edge detection, were used for corneal segmentation. After corneal segmentation, a region of interest (ROI) corresponding to central corneal stroma was selected (Fig. 1). The mean number of pixels corresponding with corneal tissue was approximately 9500 pixels, of which approximately one-third (3000 pixels) were included in the ROI. Image pixels corresponding to corneal epithelium were omitted from the ROI as they carried statistical information different from that of the stroma.²⁰ The width of the ROI was selected for each participant, according to this principle. The length of ROI (200 pixels) was optimized to achieve the highest range of changes in the statistical parameters of corneal deformation during the air-puff measurement; that is, to maximize the amount of signal-carrying data. Despite the relatively small variation in ROI width, the ROI in each image was normalized by dividing by the number of non-zero pixels, to avoid potential bias due to differences in the number of pixels within images. 

Figure 1

 

Corneal Scheimpflug images and the ROI (yellow lines) selected for analysis for a randomly chosen participant. (A) The shape of the cornea prior before the air-puff; (B) the shape of the cornea at the highest concavity. On the right is the corresponding histogram of pixel intensities in each ROI without normalization.

Pixels corresponding to a given ROI were grouped in a histogram according to pixel intensity (Fig. 1). The resulting histogram was further approximated by the Gamma probability density function,  

\(\def\upalpha{\unicode[Times]{x3B1}}\)\(\def\upbeta{\unicode[Times]{x3B2}}\)\(\def\upgamma{\unicode[Times]{x3B3}}\)\(\def\updelta{\unicode[Times]{x3B4}}\)\(\def\upvarepsilon{\unicode[Times]{x3B5}}\)\(\def\upzeta{\unicode[Times]{x3B6}}\)\(\def\upeta{\unicode[Times]{x3B7}}\)\(\def\uptheta{\unicode[Times]{x3B8}}\)\(\def\upiota{\unicode[Times]{x3B9}}\)\(\def\upkappa{\unicode[Times]{x3BA}}\)\(\def\uplambda{\unicode[Times]{x3BB}}\)\(\def\upmu{\unicode[Times]{x3BC}}\)\(\def\upnu{\unicode[Times]{x3BD}}\)\(\def\upxi{\unicode[Times]{x3BE}}\)\(\def\upomicron{\unicode[Times]{x3BF}}\)\(\def\uppi{\unicode[Times]{x3C0}}\)\(\def\uprho{\unicode[Times]{x3C1}}\)\(\def\upsigma{\unicode[Times]{x3C3}}\)\(\def\uptau{\unicode[Times]{x3C4}}\)\(\def\upupsilon{\unicode[Times]{x3C5}}\)\(\def\upphi{\unicode[Times]{x3C6}}\)\(\def\upchi{\unicode[Times]{x3C7}}\)\(\def\uppsy{\unicode[Times]{x3C8}}\)\(\def\upomega{\unicode[Times]{x3C9}}\)\(\def\bialpha{\boldsymbol{\alpha}}\)\(\def\bibeta{\boldsymbol{\beta}}\)\(\def\bigamma{\boldsymbol{\gamma}}\)\(\def\bidelta{\boldsymbol{\delta}}\)\(\def\bivarepsilon{\boldsymbol{\varepsilon}}\)\(\def\bizeta{\boldsymbol{\zeta}}\)\(\def\bieta{\boldsymbol{\eta}}\)\(\def\bitheta{\boldsymbol{\theta}}\)\(\def\biiota{\boldsymbol{\iota}}\)\(\def\bikappa{\boldsymbol{\kappa}}\)\(\def\bilambda{\boldsymbol{\lambda}}\)\(\def\bimu{\boldsymbol{\mu}}\)\(\def\binu{\boldsymbol{\nu}}\)\(\def\bixi{\boldsymbol{\xi}}\)\(\def\biomicron{\boldsymbol{\micron}}\)\(\def\bipi{\boldsymbol{\pi}}\)\(\def\birho{\boldsymbol{\rho}}\)\(\def\bisigma{\boldsymbol{\sigma}}\)\(\def\bitau{\boldsymbol{\tau}}\)\(\def\biupsilon{\boldsymbol{\upsilon}}\)\(\def\biphi{\boldsymbol{\phi}}\)\(\def\bichi{\boldsymbol{\chi}}\)\(\def\bipsy{\boldsymbol{\psy}}\)\(\def\biomega{\boldsymbol{\omega}}\)\(\def\bupalpha{\unicode[Times]{x1D6C2}}\)\(\def\bupbeta{\unicode[Times]{x1D6C3}}\)\(\def\bupgamma{\unicode[Times]{x1D6C4}}\)\(\def\bupdelta{\unicode[Times]{x1D6C5}}\)\(\def\bupepsilon{\unicode[Times]{x1D6C6}}\)\(\def\bupvarepsilon{\unicode[Times]{x1D6DC}}\)\(\def\bupzeta{\unicode[Times]{x1D6C7}}\)\(\def\bupeta{\unicode[Times]{x1D6C8}}\)\(\def\buptheta{\unicode[Times]{x1D6C9}}\)\(\def\bupiota{\unicode[Times]{x1D6CA}}\)\(\def\bupkappa{\unicode[Times]{x1D6CB}}\)\(\def\buplambda{\unicode[Times]{x1D6CC}}\)\(\def\bupmu{\unicode[Times]{x1D6CD}}\)\(\def\bupnu{\unicode[Times]{x1D6CE}}\)\(\def\bupxi{\unicode[Times]{x1D6CF}}\)\(\def\bupomicron{\unicode[Times]{x1D6D0}}\)\(\def\buppi{\unicode[Times]{x1D6D1}}\)\(\def\buprho{\unicode[Times]{x1D6D2}}\)\(\def\bupsigma{\unicode[Times]{x1D6D4}}\)\(\def\buptau{\unicode[Times]{x1D6D5}}\)\(\def\bupupsilon{\unicode[Times]{x1D6D6}}\)\(\def\bupphi{\unicode[Times]{x1D6D7}}\)\(\def\bupchi{\unicode[Times]{x1D6D8}}\)\(\def\buppsy{\unicode[Times]{x1D6D9}}\)\(\def\bupomega{\unicode[Times]{x1D6DA}}\)\(\def\bupvartheta{\unicode[Times]{x1D6DD}}\)\(\def\bGamma{\bf{\Gamma}}\)\(\def\bDelta{\bf{\Delta}}\)\(\def\bTheta{\bf{\Theta}}\)\(\def\bLambda{\bf{\Lambda}}\)\(\def\bXi{\bf{\Xi}}\)\(\def\bPi{\bf{\Pi}}\)\(\def\bSigma{\bf{\Sigma}}\)\(\def\bUpsilon{\bf{\Upsilon}}\)\(\def\bPhi{\bf{\Phi}}\)\(\def\bPsi{\bf{\Psi}}\)\(\def\bOmega{\bf{\Omega}}\)\(\def\iGamma{\unicode[Times]{x1D6E4}}\)\(\def\iDelta{\unicode[Times]{x1D6E5}}\)\(\def\iTheta{\unicode[Times]{x1D6E9}}\)\(\def\iLambda{\unicode[Times]{x1D6EC}}\)\(\def\iXi{\unicode[Times]{x1D6EF}}\)\(\def\iPi{\unicode[Times]{x1D6F1}}\)\(\def\iSigma{\unicode[Times]{x1D6F4}}\)\(\def\iUpsilon{\unicode[Times]{x1D6F6}}\)\(\def\iPhi{\unicode[Times]{x1D6F7}}\)\(\def\iPsi{\unicode[Times]{x1D6F9}}\)\(\def\iOmega{\unicode[Times]{x1D6FA}}\)\(\def\biGamma{\unicode[Times]{x1D71E}}\)\(\def\biDelta{\unicode[Times]{x1D71F}}\)\(\def\biTheta{\unicode[Times]{x1D723}}\)\(\def\biLambda{\unicode[Times]{x1D726}}\)\(\def\biXi{\unicode[Times]{x1D729}}\)\(\def\biPi{\unicode[Times]{x1D72B}}\)\(\def\biSigma{\unicode[Times]{x1D72E}}\)\(\def\biUpsilon{\unicode[Times]{x1D730}}\)\(\def\biPhi{\unicode[Times]{x1D731}}\)\(\def\biPsi{\unicode[Times]{x1D733}}\)\(\def\biOmega{\unicode[Times]{x1D734}}\)\begin{equation}f\left( x \right) = {{{\beta ^\alpha }} \over {{\rm{\Gamma }}\left( \alpha \right)}}{x^{\alpha - 1}}{e^{ - \beta x}},\end{equation}

where α is the shape parameter, β is the rate parameter, and Display Formula\({\rm{\Gamma }}\left( \alpha \right)\) is the standard Gamma function. These parameters were estimated using the method of maximum likelihood from the pixel intensities of the selected ROI in each image, in a way similar to that of the Generalized Gamma distribution performed previously for static OCT images.²⁰ In this way, the two parameters form two time series for the duration of each measurement. The actual optimization approach applied for fitting intensity values to the Gamma distribution was done using the ‘gamfit' build-in Matlab function (version 9.4 [R2018a]; MathWorks, Natick, MA, USA). The procedure was as follows. First, the method of moments was used as a statistical method for estimation of the initial starting point for the maximum likelihood estimates. Second, it was ensured that maximum likelihood was possible, otherwise the first step had to be repeated. Finally, the maximum likelihood estimates were found.²⁵ The basic steps of the methodology to obtain α and β parameters from each frame are shown in Figure 2.  

Figure 2

 

Main steps to obtain α and β parameters from each frame. Step 1: corneal segmentation from the original corneal image, the blue and purple lines represents the anterior and posterior surfaces of the cornea, respectively. Step 2: ROI selection (a region of corneal stroma). Step 3: Histogram that represents the intensity of the pixels of ROI. Step 4: Probability density function (PDF) of the Gamma function, represented by the red line, fitted to the pixels in ROI. The fit is achieved by estimating the two parameters of the Gamma distribution (α and β), estimated using the method of maximum likelihood.

After estimation of α and β parameters was complete, an eighth order polynomial was fit to the time series of the parameter α(t) to find the inflexion points and use them to distinguish five phases of corneal deformation during the measurement:² cornea at rest 1, that is, cornea in the convex state before air-puff stimulation (I_a); applanation 1, that is, first transition zone (II_a); concave cornea (III); applanation 2, that is, second transition zone (II_b); and cornea at rest 2, that is, cornea in convex state after air-puff stimulation (I_b), as is shown in Figure 3. The time series of the α(t) and β(t) parameters was further analyzed using 10 parameters (P1, P2,…, P10) described in Table 1 and further included for statistical analysis. 

Figure 3

 

Example of the dynamic evolution of the shape (α) and rate (β) parameters extracted from the fitting the Gamma distribution to the light intensity of the cornea over 140 images (frames) acquired for a single Corvis ST measurement for one measurement. The orange dashed line shows the fit of an eighth order polynomial to the α(t) parameter. This procedure is necessary to demarcate the five areas: cornea in convex state before air-puff stimulation (I_a), first transition zone (II_a), cornea in concave state (III), second transition zone (II_b), and cornea in convex state after air-puff stimulation (I_b).

Table 1

 

Description of the Parameters Related to Corneal Microstructure

 

Sstatistical analysis was performed using Microsoft Office Excel (Microsoft Office Professional Plus 2013; Microsoft; Redmond, WA, USA). The Shapiro-Wilk test was used to test the distribution type (Gaussian or non-Gaussian) of all continuous variables. The independent two-sample t-test was used to assess differences in macroscopic and microscopic biomechanical parameters between the KC and control groups. When the condition of normality was not met, the Mann-Whitney U test was performed instead. The level of significance was set to 0.05 and a Bonferroni correction was applied to address the problem of multiple comparisons. 

Results

The microscopic biomechanical parameters α(t) and β(t) were good discriminators between keratoconus and control eyes (Fig. 4), with statistically significant differences for both α(t) (Mann-Whitney U test, P < 0.05/N [N = 140, Bonferroni]) and β(t) (Mann-Whitney U test, P < 0.05/N [N = 140, Bonferroni]). Similarly, the derived parameters from the time evolution of α(t) and β(t), (P1,…, P10), described in Table 1, also were good discriminators between keratoconus and control eyes (Table 2). 

Figure 4

 

Mean values of the dynamic evolution of the shape, α(t), and rate, β(t), parameters extracted from the fitting of Gamma distribution to the pixel intensities of corneal Scheimpflug images for 10 keratoconus participants (blue) and 10 control participants (red). Error bars: SE at 95% confidence level.

Table 2

 

Comparison of Biomechanical Macroscopic Parameters Extracted From the Corvis ST Software^17,27 and Biomechanical Microscopic Parameters, Defined in Table 1, Extracted From Statistical Modeling of Corneal Light Intensity Distribution, Between the KC (N = 10) and Control (N = 10) Groups

 

The macroscopic and microscopic biomechanical parameters with the best discriminating potential between keratoconus and control eyes were extracted from Table 2, defined as those parameters with the smallest overlapping ranges between both groups. Using that criterion, the best macroscopic parameters were the maximum deformation amplitude, central corneal thickness (CCT), and IOP, as acquired with the Corvis ST, and the best microscopic parameters were P2, P3, and P4 (Table 1). Plotting combinations of macroscopic and microscopic parameters formed the most successful discriminating tool (Figs. 5G–I), which performed better than plots of only microscopic parameters (Figs. 5D–F) or only macroscopic parameters (Figs. 5A–C). 

Figure 5

 

Scatter plots demonstrating the discriminative power between keratoconus (N = 10) and control (N = 10) eyes of: (A–C) traditional macroscopic biomechanical parameters extracted from Corvis ST software, that is, maximum deformation amplitude, CCT, and IOP; (D–F) novel microscopic biomechanical parameters extracted from the proposed image processing methodology, the definition of parameters P2, P3, P4 is described in Table 1; (G–I) combination of the previously shown macroscopic and microscopic biomechanical parameters. Gray dashed lines in (G–I) indicate the separation between keratoconus and control eyes.

Discussion

This study derived a set of novel parameters from the light intensity distribution of Scheimpflug images of the corneal stroma during mechanical stimulation. These parameters relate to the microstructure of corneal tissue^20,24 and show great potential to discriminate between keratoconus and healthy eyes (Fig. 4). Nevertheless, a combination of traditional macroscopic biomechanical parameters and the newly introduced microscopic biomechanical parameters was the most successful in discriminating between keratoconus and healthy controls (Figs. 5G–I). 

Different clinical studies attempted to use macroscopic biomechanical parameters taken from Corvis software as a discriminating tool,^14,16,17,26 but the results often were inconclusive, especially in early keratoconus cases.^16,26 The results of our macroscopic biomechanical parameters in Table 2 are in agreement with those of previous clinical works.^14,16,17,26 Moreover, a combination of macroscopic Corvis parameters was proposed by Vinciguerra et al.²⁷ and added to device software under the name of the Corvis Biomechanical Index (CBI) as a risk assessment tool for ectasia. CBI ranges from 0.00 (low risk of ectasia) to 1.00 (high risk of ectasia). Although all keratoconic cases in this study were correctly classified by the CBI as keratoconus, 3 of 10 controls scored over 0.80, despite no signs of ectasia in the clinical examination. These false-positives indicated the limitations of a diagnosis based solely on macroscopic parameters. 

Combining macroscopic and microscopic parameters provided the best level of discrimination between keratoconus and controls (Figs. 5G–I). Furthermore, microscopic parameters extracted from image processing on their own were better discriminators than the macroscopic parameters on their own (Figs. 5A–F). Regardless of the limitations of the geometric parameters extracted from corneal dynamic stimulations, the deformation amplitude was repeatedly the best isolated discriminant parameter.^14,16 Similarly, CCT is considered by ophthalmologists as a key parameter in keratoconus diagnosis. However, it is noteworthy that the combination of these two macroscopic parameters (Fig. 5A) leads to some misclassifications, while the combination of these macroscopic parameters with microscopic parameters obtained from image statistical analysis did not (Figs. 5G, 5H). In addition, the differences found between normal and keratoconus eyes before mechanical excitation of the cornea (the first frames in Fig. 4) suggested that the observed differences were related to inner microstructural differences between keratoconus and normal cornea rather than just being an artefact created by light scattering during corneal mechanical stimulation. The differences in α and β parameters between keratoconic and control eyes at the very first frame (t₀ = 0 s) were statistically significant (t-test, for α(t₀): P = 0.004; for β(t₀): P = 0.003). In addition this suggested that the image processing method presented also could be applied to static Scheimpflug images. Therefore, the α and β parameters could assist practitioners with keratoconus diagnosis based on Scheimpflug topography, which is amply used in clinical settings. Corneal densitometry, defined as a map of the amount of backscattered light in the different regions of the cornea,²⁸ is provided as an add-on to the standard software of commercially available Scheimpflug tomography. We would expect our results to be correlated with traditional densitometry values, since both methods are based on light backscattering, but this is yet to be confirmed. 

A general interpretation of the physical meaning of α and β parameters when fitting Gamma function to the light intensity distribution of pixels is available in the literature.^20,29 In particular, it has been suggested that the scale parameter α quantifies the average backscatter power and can be related to the average scattered cross-section, whereas the shape parameter β was related to the scatter density. This study does not provide sufficient data to assess the origin of the observed differences in pixel intensity distributions between control and keratoconus eyes. However, based on previous studies^20,29 and the well-known fact that in keratoconus corneal biomechanics is compromised,^11–13 it is possible to make an educated guess that the underlying mechanism is based on differences in distribution of the tissue scatter due to microstructural alterations caused by the disease process. 

The small sample size could be seen as a limitation of this proof-of-concept study. However, it is noteworthy that statistically significant differences were found despite this fact or the relatively small number of pixels in each image (200 × 576 pixels compared to, for example, 600 × 800 pixels in the Pentacam HR). Also, it is noteworthy that the keratoconic cases included in the study were rather mild (as moderate or advanced cases were excluded), strengthening the potential of the presented results. 

In conclusion, to our knowledge, we present the first set of parameters related to the corneal microstructure, derived from statistical modeling of the light intensity distributions of dynamic Scheimpflug images. This novel approach showed the potential of combining macroscopic and microscopic corneal properties to successfully discriminate between mild keratoconus and control eyes. 

Acknowledgments

Supported by Grant 779960 from the European Union's Horizon 2020 research and innovation program, the Statutory Funds of Wroclaw University of Science and Technology, and Grant 666295 (MW) from the European Union's Horizon 2020 research and innovation program. 

Disclosure: A. Consejo, None; K. Gławdecka, None; K. Karnowski, None; J. Solarski, None; J.J. Rozema, None; M. Wojtkowski, None; D.R. Iskander, None