Two-Dimensional Plane for Multi-Scale Quantification of Corneal Subbasal Nerve Tortuosity

Abstract

Purpose: To assess the performance of a novel system for automated tortuosity estimation and interpretation.

Methods: A supervised strategy (driven by observers' grading) was employed to automatically identify the combination of tortuosity measures (i.e., tortuosity representation) leading to the best agreement with the observers. We investigated 18 tortuosity measures including curvature and density of inflection points, computed at multiple spatial scales. To leverage tortuosity interpretation, we propose the tortuosity plane (TP) onto which each image is mapped. Experiments were carried out on 140 images of subbasal nerve plexus of the central cornea, covering four levels of tortuosity. Three experienced observers graded each image independently.

Results: The best tortuosity representation was the combination of mean curvature at spatial scales 2 and 5. These tortuosity measures were the axes of the proposed TP (interpretation). The system for tortuosity estimation revealed strong agreement with the observers on a global and per-level basis. The agreement with each observer (Spearman's correlation) was statistically significant (α_s = 0.05, P < 0.0001) and higher than that of at least one of the other observers in two out of three cases (ρ_OUR = 0.7594 versus ρ_Obs3 = 0.7225; ρ_OUR = 0.8880 versus ρ_Obs1 = 0.8017, ρ_Obs3 = 0.7315). Based on paired-sample t-tests, these improvements were significant (P < 0.001).

Conclusions: Our automated system stratifies images by four tortuosity levels (discrete scale) matching or exceeding the accuracy of experienced observers. Of importance, the TP allows the assessment of tortuosity on a two-dimensional continuous scale, thus leading to a finer discrimination among images.

Changes in corneal nerves have been observed in various conditions such as ocular allergies and dry eye disease, neurotrophic keratopathy, infectious keratitis, contact lens wear, keratorefractive surgeries, corneal transplantation, and cataract surgery.^1–7 For clinical evaluation of the corneal nerve structures in patients, in vivo confocal microscopy (IVCM) is often used, which visualizes corneal nerves and in particular the subbasal nerve plexus.^6,8,9 In addition to measuring the corneal nerve density, this imaging modality enables the assessment of morphologic features of the nerves such as width, reflectivity, orientation, branching patterns, bead-like formations and tortuosity.^10,11 

Increased tortuosity, in particular, has been associated with various conditions, such as dry eye disease, Graves' ophthalmopathy, corneal infections, small-fiber neuropathy, and diabetic retinopathy.^12–19 However, tortuosity has no widely accepted definition, and it is often estimated qualitatively by comparing patient images with reference images.¹¹ Moreover, several corneal nerve fibers with different tortuosity characteristics (e.g., amplitude, number of inflection points) might be present in the same image. This inevitably introduces a degree of subjectivity; therefore, quantitative, objective and repeatable methods of tortuosity estimation would advance the field. 

Several quantitative measures of tortuosity have been proposed for the retinal vasculature as observed in fundus camera images,^20–27 and for the brain vasculature²⁸ (see Ref. 29 for a recent comparative study). In contrast, there is very limited published work on corneal nerve fibers; Kallinikos et al.³⁰ were the first to propose an objective, semiautomated method for quantifying the tortuosity of corneal nerve fibers. Scarpa et al.³¹ adapted the algorithm proposed by Grisan et al.²¹ for retinal vessels to work with corneal nerve fibers. 

Previously reported approaches have the following limitations, which we address herein: (1) they are single-scale, as they do not differentiate between high- and low-frequency turns often characterizing corneal nerve fibers (Fig. 1); (2) they are based on a single tortuosity index (or coefficient) in which the combination and weight of each tortuosity measure (e.g., curvature, number of inflection points) is defined a priori; (3) the estimated tortuosity consists of a single number, therefore they are lacking in interpretation as the contribution of each tortuosity measure in the final estimate is lost; (4) they do not provide the confidence level of the estimated tortuosity, thus ignoring that some tortuosity levels might be more difficult to distinguish from others. 

Figure 1

 

Example of multiple spatial scales observed in a corneal nerve fiber (a, yellow arrows). Multi-scale decomposition (c, d) of a synthetic corneal nerve fiber (b).

In our preliminary work on corneal nerve fiber tortuosity estimation,³² we addressed points 1 and 2, investigating which role is played by different spatial scales (i.e., different frequency of turns along the fiber) in tortuosity estimation (Fig. 1). Experimental results suggested that the notion of tortuosity can be captured quantitatively in a richer and more informative way by computing the tortuosity measures at multiple spatial scales, and we proposed an automated methodology to compute it. The need to measure tortuosity at multiple spatial scales is further strengthened by a recent study³³ highlighting the roles played by “short- and long-range tortuosity” in tortuosity grading, which we (a) formalized mathematically with multi-scale analysis and (b) discovered through a machine learning algorithm.³² 

Herein, we assess the performance of our automated methodology in greater depth and compare it with the tortuosity assessments provided by expert observers. We measure performance on a global and per-level basis, to evaluate the accuracy of our system in predicting correctly all the tortuosity levels and each tortuosity level, respectively. Moreover, we address points 3 and 4 by proposing the tortuosity plane (TP), a two-dimensional (2-D) plane in which the coordinate axes are the weighted mean curvatures at the most significant spatial scales identified automatically. In this framework, each IVCM image is represented as a point on the TP, thus allowing a better interpretation of the estimated tortuosity. Of importance, the TP represents images on a continuous 2-D tortuosity scale, thus leading to a finer discrimination among images. 

Methods

Image Acquisition

One hundred forty de-identified images of the subbasal nerve plexus of the central cornea were selected from an existing database. Images were acquired from 76 subjects (49 women and 27 men, 54.9 ± 19.3 years old) with laser scanning IVCM (Heidelberg Retina Tomograph 3 with the Rostock Cornea Module; Heidelberg Engineering GmbH, Heidelberg, Germany). The diode laser has a 670-nm red wavelength, and the microscope is equipped with a ×63 objective lens with a numerical aperture of 0.9 (Olympus, Tokyo, Japan). The images obtained represent a coronal section of the cornea of 400 × 400 μm², from the subbasal layer. For imaging, a disposable sterile polymethylmethacrylate cap (Tomo-Cap; Heidelberg Engineering) was filled with a layer of hypromellose 0.3% gel (GenTeal; Alcon, Fort Worth, TX, USA) and mounted on the Cornea Module. After topical anesthesia with 0.5% proparacaine hydrochloride (Alcon), a drop of hypromellose 0.3% gel was applied to both eyes. One drop of this gel was also put on the outside tip of the cap to improve optical coupling. Then, the Cornea Module of the confocal microscope was advanced until the gel contacted the central surface of the cornea. Digital images were then recorded at a rate of 3 frames per second in sequence mode, including 100 images per sequence, from all layers of the cornea. 

For this pilot-level study, 140 images showing a wide range of tortuosity characteristics were chosen from patients with a variety of eye conditions, typically at a depth of 50 to 80 μm (Fig. 2). Specifically, 20 images were taken from healthy subjects, 100 images from patients with dry eye disease, and 20 images from subjects with herpes simplex keratitis. Images were graded independently by three expert observers (PH, AK, and SA) using four tortuosity levels according to a previously published protocol.¹¹ Nerves were traced using NeuronJ (in the public domain at http://www.imagescience.org/meijering/software/neuronj/), a plug-in for the ImageJ software (in the public domain at http://imagej.nih.gov/ij). 

Figure 2

 

Corneal nerve fibers in the subbasal layer in IVCM. Original images from tortuosity level 1 (A) to level 4 (D).

The study was performed in accordance with the tenets of the Declaration of Helsinki. 

Proposed Tortuosity Estimation and Interpretation System

Our methodology can be summarized in four steps: (1) to compute multiple-scale tortuosity measures (or features) from fiber centerlines; (2) to use a feature selection approach to identify automatically the combination of tortuosity measures and scales leading to the best agreement with experienced observers; (3) to assign a new image to one of the four tortuosity levels according to the value of its tortuosity measures; and (4) to map the new image as a point on the TP and show the level of confidence (i.e., reliability) for the tortuosity level assigned by the system. 

MATLAB (Release 2014b, The MathWorks, Inc., Natick, MA, USA) was used to implement the proposed system. 

Computing Multiple-Scale Shape Features

The initial feature pool is formed by 18 measurements taken on the fibers (Table 1). We use maximum curvature and inflection points, the basis of the best-performing tortuosity indices in recent reports (see introduction), and average curvature, which led to promising performance in several methods.^22,24,25,34 These three measures are computed at six spatial scales to take into account the contribution of turns at six different frequencies. As Figure 3 shows, we start by measuring curvatures and inflections points on the original fiber (including all the frequencies, scale 1). Then, we smooth out its shape to remove high-frequency turns and concentrate on lower frequency turns (scales 2–6). 

Table 1

 

Tortuosity Measures Investigated With Our Automated Feature Selection Approach*

 

Figure 3

 

Multi-scale analysis of a corneal nerve fiber. From left to right: the original fiber including turns at all frequencies (spatial scale 1), and its smoothed versions at spatial scales from 2 to 6 to take into account turns at intermediate and high frequencies.

Curvature and inflection points are measured with a highly accurate algorithm specifically designed to work on anatomic structures imaged at low resolution (such as corneal nerve fibers in our data set).³² Notice, we use the density of inflection points (i.e., number of inflection points per unit length) rather than their number, to yield a consistent comparison among fibers of different length. 

Since each IVCM image includes a variable number of corneal nerve fibers, we average the 18 measures over all fibers in the image to obtain 18 global (i.e., image-level) measures. Assuming that longer fibers are more informative than shorter ones, we apply weighted averaging in which each fiber contribution to the global measure is weighted by its length. 

Identifying the Best Combination of Tortuosity Measures

Feature selection^35–37 is a class of machine learning algorithms identifying the most discriminative features (or measures) for a task, here for tortuosity assessment. As our pool of measures is computed at multiple spatial scales, it also selects the most discriminative scales. 

Here, we use a wrapper strategy,³⁷ in which the classification is based on a multinomial logistic ordinal regressor (MLOR). The MLOR models the log-cumulative odds³⁸; that is, the logarithm of the ratio of two probabilities: the probability that an image y has level of tortuosity c_j or lower, Pr(y ≤ c_j) and the probability that its tortuosity level is greater than c_j, Pr(y > c_j), using a linear combination of the p tortuosity measures (features) X_i:  

 

where, j ∈ {1,2,3} and α and β are the weights of each tortuosity measure into the final assignment.  

An MLOR model is built for each combination of tortuosity measures (an optimization step is performed to estimate the weights α and β), and it is used to assign the tortuosity level to a subset of training images. In this framework, the best combination of tortuosity measures and their scale is the one minimizing the number of misclassified images. Notice, we used cross-validation (i.e., randomly splitting the data set in training and test sets) to provide a fair assessment. 

Assigning New Images to a Tortuosity Level

Once the best combination of tortuosity measures is found, we build the final MLOR model (i.e., using all the training images). The best weights α and β are estimated for each case (i.e., level 1 versus level 2, 3, or 4; level 1 or 2 versus level 3 or 4; level 1, 2, or 3 versus level 4). Then, the selected tortuosity measures are computed on the new image, and its tortuosity level is obtained using the maximum a posteriori criterion: the tortuosity level that the MLOR predicts as the most likely among the four considered is assumed as the estimated tortuosity level. 

Tortuosity Plane and Confidence of the Estimated Tortuosity Level

Previously proposed methods typically provide a single tortuosity level (or coefficient value), obfuscating the different effects of the factors considered, and limiting interpretations that could be important to the ophthalmologist for diagnosis. Using the geometric interpretation of the MLOR model, we map each IVCM image onto a plane whose axes are the best tortuosity measures identified automatically by the feature selection procedure (in our case, weighted mean curvatures at spatial scales 2 and 5 corresponding to the fourth and 13th tortuosity measures in Table 1, respectively). Geometrically, the best weights α and β of the MLOR model define the best (in the MLOR sense) linear decision boundaries separating the TP into four regions corresponding to the four tortuosity levels (black lines in Fig. 4). Any new corneal nerve image can be plotted as a point on the TP by computing the values of the two best tortuosity measures (Fig. 4 shows four examples). The region containing the new point gives the estimated tortuosity level for the corneal nerve image. 

Figure 4

 

Tortuosity estimated for four corneal nerve images by the proposed method. Images are projected onto the TP as points (markers indicate the tortuosity level assigned by the expert observers), whose coordinates are the estimated mean curvature at scale 2 (high-frequency turns) and scale 5 (low-frequency turns). The level of confidence for each tortuosity estimate is encoded with color (red and blue mean high and low confidence, respectively). Intuitively, the confidence is related to the distance from the decision boundaries (indicated in black) separating each tortuosity region.

Of importance, the TP also provides a level of confidence for the estimated tortuosity level, quantifying the reliability of the system. This confidence is given by the probability of belonging to one of the four regions (i.e., tortuosity levels) estimated automatically by the MLOR model and is intuitively proportional to the distance of the point (i.e., image) from the linear decision boundaries. In fact, the closer the point to the boundary between two adjacent regions, the less reliable is the estimated tortuosity level. The level of confidence for all the points on the TP can be estimated once the system is trained (i.e., before analyzing the target images) and color coded for immediate, intuitive visualization (Fig. 4). 

Performance Criteria and Statistical Analysis

Discriminating between four tortuosity levels is a multi-class classification problem. Following the standard performance assessment methodology for these problems,³⁹ we use accuracy (Acc), sensitivity (Se), specificity (Sp), positive predicted value (Ppv), and negative predictive value (Npv), as defined below.  

 

 

 

 

 

 

 

 

  ,

where TP_i indicates true positives; TN_i, true negatives; FP_i, false positives; and FN_i, false negatives for the i-th tortuosity level (N = 4). We use weighted averaging as our data set is unbalanced, as suggested by Sokolova and Lapalme.³⁹ The above performance measures do not take into account the distance from the actual value when considering errors. To this aim, we also include the mean squared error (MSE) defined as  

 

where N_imgs is the total number of test images, y_i is the predicted tortuosity level, and l_i the true one, for the i-th image.  

In addition to the aforementioned performance measures, we use Spearman's correlation as done elsewhere (see Ref. 27 for review). To assess the statistical significance of improvements among tortuosity estimation methods we use pair-sample t-tests. 

Results

System Agreement With Consensus Ground Truth

Assessing tortuosity is a difficult and subjective task. As such, the level of agreement among expert observers is arguably smaller compared with other tasks such as density estimation (see, for example, the repeatability study by Petropoulos et al.⁴⁰) To address this problem, we create a consensus ground truth (CGT): the tortuosity level of each image is given by the majority voting (i.e., the CGT tortuosity level for an image is the one assigned independently by at least two expert observers). Consistent with the difficulty of the task, 11 images were discarded owing to complete disagreement (all assigned levels different). We trained the system using the CGT and achieved promising global-level performance, reported previously (Annunziata et al.,³² Table 1). However, some tortuosity levels might be easier to estimate compared with others, or a specific level could be predicted with poor performance. To better investigate this point, we compute each performance measure on a per-level basis and show the results in Table 2. Our system achieves almost 90% accuracy for the tortuosity levels 1 and 4. Performance decreases for the middle levels owing to smaller differences with neighboring regions. Quantitatively, sensitivity and positive-predictive values decrease considerably owing to the increase in false negatives. Qualitative considerations on the TP are made below. 

Table 2

 

Per-Level Performance Measures of the Proposed Tortuosity Estimation Approach When the CGT Is Taken as Reference

 

System Agreement With Individual Expert Observers

To explore agreement beyond consensus, we take the assessment made by one expert observer as ground truth and compare the performance of our framework (MLOR) with the other two expert observers using Spearman's correlation. This procedure is repeated three times, taking each expert observer as reference, in turn (Table 3). The correlations of our method with expert observers are significant (α_s = 0.05, P < 0.0001). Moreover, the first and second rows of Table 3 show that the correlation of our system with the reference observer is higher than that of at least one of the other observers. Computed P values for the pair-sample t-tests on the MSE show that these improvements are significant for both hypotheses “MSE_ours < MSE_Obs3” (row 1, P < 0.001) and “MSE_ours < MSE_Obs1” (row 2, P < 0.001). 

Table 3

 

Spearman's Correlation and P Values Computed When an Individual Observer (Reported as Obs1, Obs2, and Obs3) Is Taken as Reference or Ground Truth (GT)*

 

Tortuosity Plane

Figure 4 shows the TP obtained with the proposed method applied to a subset of images during a random 20-fold cross-validation procedure. The four regions related to the four tortuosity levels are identified automatically, and the best separating lines (visualized in black) are traced. Then, the two best tortuosity measures are computed on each of the testing images, and they are mapped onto the TP. Since each image belongs to a different tortuosity level, they are (correctly) mapped to four different regions. Visualizing the images as points on the TP allows an immediate and consistent interpretation. First, image 1 is close to the boundary 1-2, so the confidence of tortuosity level 1 (instead of 2) is limited; this is consistent, as some fibers appear rather tortuous for level 1. The same considerations apply to image 3. Second, although image 2 belongs to the tortuosity level with the lowest overall accuracy (level 2), it is mapped to the center of the corresponding TP region, and therefore, the confidence of the prediction is relatively high. Image 4 is classified with the highest confidence as it clearly shows an abrupt large change (scale 5) in the direction of the second fiber (counting from left to right), and also high-frequency changes in other fibers (scale 2). Previous, single-scale methods would not tease out this difference. Third, the contribution of high- and low-frequency turns (i.e., curvatures at scales 2 and 5) for the predicted tortuosity level can be easily deduced. For instance, images 1 and 2 show a relatively low mean curvature at scale 5, but this is higher for images 3 and 4. Then, although image 2 is correctly assigned to tortuosity level 2, a higher mean curvature at scale 2 than image 3 is present: this suggests that, in our data set, the mean curvature at scale 2 (i.e., high-frequency turns) is slightly less discriminative than the one at scale 5. 

We ran a leave-one-out cross-validation simulating the scenario in which the tortuosity of a new unseen image is assessed. For each validation, the best boundaries were estimated automatically, and the tortuosity level was assigned. Figure 5 shows the result of this process. We use different markers to indicate the true tortuosity level (i.e., l in Equation 7) of each image as per CGT. The four clusters reflecting the tortuosity levels are clearly visible and separated automatically by our system, although some misclassifications are present. We notice that the maximum error made by the system is always within one tortuosity level, showing very low MSE errors. Moreover, different regions show different spreading on the TP, especially in region 2 (very low) and region 4 (very high), suggesting that intra-level tortuosity characteristics should be better investigated in future studies. Of interest, tortuosity does not vary much in the Y axis, but does vary in the X axis, for levels 1 and 2. At levels 3 and 4, tortuosity varies in both the Y and X axes. 

Figure 5

 

All IVCM images are projected onto the TP after a leave-one-out cross-validation on unseen images. Markers indicate the CGT tortuosity level. The estimated tortuosity level corresponds to the region within which an IVCM image is mapped. During each cross-validation, the best decision boundaries (shown in black) are computed based on the training set.

Discussion

We have presented a novel methodology, the tortuosity plane, to represent and quantify tortuosity in the context of corneal nerve fibers in IVCM images of the subbasal nerve plexus. Unlike other methods using a single index in which the combination and weight of each tortuosity measure are defined a priori, we capture tortuosity with two measures selected by an automatic procedure on the basis of their agreement with grading from experienced observers. 

Our approach reveals that tortuosity is best represented with curvature-related measures taken at multiple spatial scales (specifically, scales 2 and 5 accounting for high- and low-frequency turns), a computational finding supported by a recent, general study on the perception of tortuosity of corneal nerve fibers by expert observers.³³ Unlike this previous study, the proposed method attempts to discover the best tortuosity representation and the most discriminative spatial scales for a specific data set, using a fully automated procedure (feature selection). This novel approach to tortuosity estimation has the potential to answer interesting research questions related to the role played by each tortuosity measure and spatial scale for specific disease conditions. 

Visualizing the predictions made by the proposed tortuosity estimation method on a 2-D plane, the TP, offers various advantages compared with traditional approaches providing a single number (i.e., tortuosity level or index). 

First, it allows a finer (i.e., continuous) tortuosity scale compared with methods estimating tortuosity using a few levels. Compared with approaches based on indices (i.e., providing a continuous number), the TP provides a 2-D continuous tortuosity scale, thus allowing a better interpretation of the estimated tortuosity. Moreover, the contribution of each tortuosity measure into the final estimate becomes trivial (i.e., it is readily available on the axes of the plane). Overall, the TP allows a better tortuosity stratification and the possibility to identify subcategories clustering in specific regions of the plane. 

Second, it visualizes the level of confidence (i.e., reliability) of the estimated tortuosity simply and intuitively, using color and also by the distance of a point from its closest boundaries (Fig. 4), an aspect ignored so far. As Figure 4 shows, the proposed method takes into account that the tortuosity of some images can be estimated with less confidence compared with others. A visual investigation suggests that a decreased level of confidence is due to the presence of fibers with different tortuosity characteristics within the same image. This is the result of training our system with image-level information (expert observers were asked to provide the level of tortuosity for the whole image). As expected, tortuosity levels 1 and 4 are estimated with higher confidence on average (e.g., regions 1 and 4 are mainly colored in red), compared with levels 2 and 3. 

We showed in a pilot experiment that images with four tortuosity levels are grouped in distinct clusters on the TP. This immediately suggests that the TP could provide a quantitative, intuitive tool to stratify patients by condition, or by stage of a specific disease, a research avenue that we are pursuing. 

Quantitative results suggest that our system achieves a level of performance that in two out of three cases is significantly better than that of experienced observers when compared against each other. Notice that we assessed the performance in a more challenging experimental setting compared with previous approaches (i.e., using four tortuosity levels instead of, at most, three in previous work). 

The aforementioned considerations demonstrate the importance of using a 2-D plane and a multi-scale tortuosity estimation over single tortuosity indices, as proposed previously. Notice that the best scale for each tortuosity measure is selected automatically by our feature selection strategy. This is important as it avoids taking into account tortuosity measures computed at very high scales (i.e., very low frequency). In fact, these tortuosity measures could be affected by the presence of branching points breaking the continuity of the main fiber. We found no significant correlation between the tortuosity assigned by our system and the number of branches per image (ρ = 0.0174, P = 0.85). Likewise, averaging each tortuosity measure at fiber level is important to avoid the effect of higher-order aberration, which could favor longer fibers, for instance. No significant correlation between total nerve length per image and the tortuosity assigned by the system was found (ρ = −0.1693, P = 0.0657). 

However, our method has some limitations: (1) Measurements are taken on manually traced fibers, thus introducing a certain level of subjectivity. Combining the method proposed herein with a reliable algorithm for automated fiber tracing would eliminate this residual subjectivity and make the whole system fully automated. (2) Since corneal nerve images were annotated as a whole (i.e., each image is associate with a single tortuosity level by the experienced observers), the contribution of each fiber is lost. We assumed that longer fibers should weight more in the global tortuosity measure, but other possibilities can be explored. (3) The number of tortuosity measures leading to the best agreement with experienced observers depends on the data set used, and more than two measures could lead to better performance on other data sets. This poses a visualization question: three tortuosity measures can still be visualized in a tortuosity volume, but four or more cannot. Although further tests are needed, the need for a higher dimensional representation including more than three tortuosity measures seems unlikely, given that we used a representative set of clinical IVCM images including different diseases. 

Our immediate future work will explore the combination of this method for automated tortuosity estimation and interpretation with automated corneal nerve fiber detection.^41–43 This will lead to a fully automated system, which would contribute to disease stratification by tortuosity and better diagnosis based on objective, highly informative measurements. 

Acknowledgments

Supported by European Union Marie Sklodowska-Curie Initial Training Network Retinal Vascular Modelling, Measurement And Diagnosis (REVAMMAD), Project Number 316990 (RA, ET) and by the National Institutes of Health Grant NIH R01-EY022695 (PH), Falk Medical Research Foundation (PH), and Massachusetts Eye and Ear Infirmary Foundation (PH). The authors alone are responsible for the content and writing of this paper. 

Disclosure: R. Annunziata, None; A. Kheirkhah, None; S. Aggarwal, None; B.M. Cavalcanti, None; P. Hamrah, None; E. Trucco, None