Automatic corneal nerve fiber segmentation and geometric biomarker quantification

Geometric and topological features of corneal nerve fibers in confocal microscopy images are important indicators for the diagnosis of common diseases such as diabetic neuropathy. Quantitative analysis of these important biomarkers requires an accurate segmentation of the nerve fiber network. Currently, most of the analysis are performed based on manual annotations of the nerve fiber segments, while a fully automatic corneal nerve fiber extraction and analysis framework is still needed. In this paper, we establish a fully convolutional network method to precisely enhance and segment corneal nerve fibers in microscopy images. Based on the segmentation results, automatic tortuosity measurement and branching detection modules are established to extract valuable geometric and topological biomarkers. The proposed segmentation method is validated on a dataset with 142 images. The experimental results show that our deep learning-based framework outperforms state-of-the-art segmentation approaches. The biomarker extraction methods are validated on two different datasets, demonstrating high effectiveness and reliability of the proposed methods.


Introduction
Corneal confocal microscopy (CCM) is an efficient and non-invasive imaging technique that is used to examine the human corneal nerve fiber morphology in a variety of diseases  3]. There exist strong connections between disease progression and geometric/topological changes of nerve fibers. Previous findings show that tortuosity, branching points and density, length and angles of corneal nerve fibers are potential indicators to reflect the function of neuropathy severity [4][5][6][7][8] and type-II diabetes [1,9].
Diabetic neuropathy is one of the common long-term complications of diabetes, in which the nerve fibers in the body are affected due to hyperglycemia [2]. It includes a set of neuropathy such as peripheral neuropathy [3] and autonomic neuropathy [10]. The diagnosis of diabetic neuropathy relies on proper measures such as nerve biopsy and skin biopsy based on the presence of symptoms such as the feeling of burning, tingling, weakness or pain in the limbs [11]. These symptoms are not always the indication of nerve damage, which makes the early detection of neuropathy very challenging. The nerve plexus innervated in the peripheral corneal is one of the most sensitive tissues of the human body providing sensors for touching, temperature and chemicals, and plays an important role in eye protection [12]. Recent studies have highlighted that impaired corneal ultrastructures are observed in patients with diabetic neuropathy [4,13,14]. Moreover, it is reported that corneal nerve loss was observed in the early development of neuropathy [9], and in some cases even before the development of other complications such as diabetic retinopathy, suggesting that corneal nerve loss might be an early indicator for diabetic neuropathy and diabetic retinopathy. The sensory nervous system in the limbs and the organs is difficult to observe and examine, while the corneal nerve fibers can be observed by CCM, which provides the potential for building a computer-aided diagnosis (CAD) system. Recently, different CAD systems [15][16][17][18][19][20][21] have been well developed in radiology. To better analyze nerve fibers, we also need to set up a regular work flow: nerve fiber segmentation, biomarker extraction and performing pathological study by investigating the association between the changes of biomarkers and the progression of the disease. However, the effectiveness of biomarker extraction relies on an accurate nerve fiber segmentation process, which is still strongly needed for the further analysis of corneal confocal microscopy images.
In general, the corneal nerve fiber segmentation can be considered as a problem within the scope of curvilinear structure enhancement and segmentation tasks [22][23][24]. For instance, extensive conventional approaches have been applied to segment elongated medical imaging structures in three categories: classifier based, tracking based and filter based. Classifierbased methods extract a set of features from image pixels and use a pre-trained classifiers to discriminate vascular pixels from background pixels. Soares et al. [25] extract a feature vector from the pixel intensity and matched filter responses to segment retinal blood vessels using a supervised classifier. Tracking-based methods iteratively grow the model from one seed point inside the structure to another seed point located at the peripheral area [26,27] to trace all vessel regions. Filter-based algorithms are based on convolution kernels which can be used to enhance curvilinear structures in the image domain. These approaches are generally faster and simpler than supervised methods. Mendonca et al. [28] employed differential filters to detect vessel centerlines followed by morphological operators for vessel segmentation. Zhang et al. [29] proposed a crossing-preserving segmentation approach, which designs multi-scale rotating filters in invertible orientation scores to highlight and extract curvilinear structures.The orientation score based framework will be further employed in this work to extract branches from corneal nerve fiber images. Al-Fahdawi et al. [30] proposed a robust and fast nerve segmentation and morphometric parameter quantification system for corneal confocal microscope images. Three main segmentation steps including anisotropic diffusion filtering based noise reduction, morphological processings and edge detection were applied to detect all the nerve fibers.
Convolutional neural networks (CNN) have gained a lot of attention in the recent years because of the effective convolution kernels which are directly learned from the training images, instead of relying on user-defined features [31][32][33][34]. One way of segmenting objects is to create a patch around each sample and apply CNN framework to classify these patches. The drawback of this scheme is that it does not make use of the larger context information around the samples. In 2015, Ronneberger et al. [35] proposed a fully convolutional network (FCN) with a U-shape (U-Net) which effectively deals with segmentation problem [36]. A U-Net-based CNN's architecture consisting of a contracting path is proposed by Colonna et al. [37] for the fully automatic tracing of corneal nerves. It demonstrates the potential of CNN in identifying clinically useful features. These deep learning-based segmentation techniques show great potential in extracting useful curvilinear structures from different medical images.
In this work, we aim to set up a complete framework to enhance and segment corneal nerve fiber structures from CCM images, and then design specific modules to extract potentially important nerve fiber biomarkers for further analysis. Therefore, in the first step, we establish a FCN-based segmentation approach to precisely extract nerve fibers and generate corresponding centerline maps. Based on the segmented nerve fiber centerline, we then designed tortuosity measurement and branch detection modules to extract never fiber biomarkers for further geometric and topological analysis. The proposed framework is evaluated on two different datasets, i.e., the corneal nerve fiber images of Maastricht (CNM) dataset and the corneal nerve tortuosity (CNT) dataset, to demonstrate the high segmentation performance and accurate biomarker measurement for further analysis.
The remainder of this paper is organized as follows. In Sect. 2, we introduce the preprocessing step for normalizing the corneal fiber images, and explain the architecture of deep learning network for the segmentation. In the second part, we show the biomarker extraction stage for tortuosity and branches. Section 3 explains and discusses the experimental results and validations for segmentation and biomarker extraction. Finally, in Sect. 5, we conclude this paper.

Methods
In this work, we first apply a normalization method to reduce the inhomogeneity in images. Then we employ a fully connected neural network architecture to segment corneal nerve fibers. Subsequently, biomarkers are extracted from the fiber segmentation results for further analysis.

Image normalization
The 2D images acquired from confocal laser scanning microscope often show significant intensity heterogeneity. The inhomogeneity can be due to photo-bleaching, attenuation in depth, image acquisition, and variations of illumination exposure rate [38]. Before we perform fiber segmentation or extract biomarkers from microscopic images, we apply image normalization to suppress intensity heterogeneity within an image and across all images.
We apply luminosity and contrast normalization developed by Foracchia et al. [39] to improve the quality of corneal nerve fibers. The method is based on the estimation of the luminosity and contrast variability in the background part of the microscopic image. The compensation of this variability is performed in the whole image. The application of this method results in remark reduction of luminosity variability and an increment of image contrast. The normalization radius for creating a disk-shaped filter is set as r LC = l h 30 , where l h is the height  Figure 1 shows one example of the effect of applying the normalization. We can observe that the inhomogeneity in the original image is suppressed effectively.

Corneal nerve fiber segmentation by U-Net
The adopted neural network in this study is based on U-Net, which had tremendous success [35] in the medical image segmentation tasks in past years. It has a U-shape and it consists of three downsampling steps (encoding path) with convolutions, followed by a batch normalization and three upsampling steps (decoding path) with repeated layers of convolution, but followed by a concatenation with the correspondingly cropped feature from contraction path and up convolutions. At the final layer, a final 1 × 1 × 1 convolution is used to map the desired number of classes, and a soft-max calculation is followed at last to output a probability for each pixel. To improve the focus on the target fiber structure of varying shapes and sizes, we incorporate attention gating (AG) modules [40] into the networks. The AG module basically performed in such a way that the input features are scaled with attention coefficient computed in AG such that the spatial regions is selected by analyzing both the activation and contextual information collected from a coarse scale. The AG module acts as memory lookup mechanisms for fine-grained information to enhance the higher level representation which is beneficial to object detection and segmentation. We used the cross entropy of a pixel-wise soft-max loss between the predicted label and ground truth label as the loss function. The architecture of the network can be found in Fig. 2.
The corneal nerve fibers in the original confocal microscope images are very thin structures with low signal-to-noise ratio and weak connections. The percentage of nerve fiber pixels only accounts for a very small portion of the whole image pixels. To obtain well-enhanced nerve fibers, a more balanced training process is designed by dilating the manually annotated nerve fiber centerlines to obtain more foreground samples in the training set. Hence, the nerve fiber maps generated by the trained model present slightly wider centerlines with better connection. In clinical study, the nerve fiber length, skeleton density, tortuosity and junction numbers are potentially more valuable biomarkers for diseases analysis compared to the width of nerve fibers. Thus, a skeletonization step is often performed after fiber extraction to obtain the centerline map of corneal nerve fibers for further biomarker analysis. This will not be affected by the increased width. To fit the data into the memory of the graphic card, every original image (with a size of 1535 × 1536) is split into smaller images with a size of 288 × 288 with overlaps. During training, data augmentation of eight times is implemented by applying a combination of random flipping and rotation of the original image considering the limited number of training images. During test, every image is also split into images with the same size and the likelihood in the overlapped area is obtained by averaging. To obtain unbiased results, we perform twofold cross validation. For a run of training, the training fold is split into a real training set (70%) and a real testing set (30%). The validation set is used to select the best classifier. We use 50 epochs for training.

Fiber tortuosity biomarker extraction
Tortuosity is a measure of the twistedness of a curve. In this study, tortuosity of corneal nerve fibers is determined using a previously described retinal vessel tortuosity index (VTI). VTI combines local and global geometric features of a vessel to quantify its tortuosity, and hence is sensitive to small tortuosity alterations. This measure of tortuosity is invariant to rigid transformations and was shown to provide a good match with visual perception of tortuosity. Mathematical derivation of VTI is given by where SD θ is the standard deviation of angles between lines tangent to the centerline at each pixel along the centerline and the x-axis. N is the number of critical points along the centerline where the first derivative of the centerline vanishes. L A and L C are lengths of the centerline and its chord, respectively. The visual demonstration of VTI is shown in Fig. 3. Finally, the magnitude M of derivation of the centerline from a straight line is obtained by averaging the ratio of centerline length to its chord length between inflection points including centerline endpoints. It is written as To automatically calculate tortuosity of the corneal nerve fibers, centerline extraction and bifurcation detection are performed on each image. MATLAB built-in "bridge" function is first used to eliminate any discontinuity along the fibers. Afterward, an iterative thinning is performed to obtain the skeleton image. Small spurs that generated during the thinning process are removed by repeatedly (20 times) deleting pixels that had only one connected nerve fiber neighborhood. The value of 20 pixels is selected because this value is approximately equal to diameter of the largest nerve fiber in the image, and the length of spurs are not expected to exceed the nerve fiber diameter. Bifurcation points are detected through convolution of the skeleton image with a unity kernel of size 3 × 3, and detection of points with a value greater or equal to 3. Nerve fiber centerlines between the bifurcation points are then labeled automatically based on their length (i.e. number of connected pixels) and tortuosity of each one is calculated correspondingly. In Sect. 4.3, we will evaluate the tortuosity measurement technique on a specific corneal nerve fiber dataset.

Detection of nerve fiber branches on orientation scores
The curvilinear topology and geometry at corneal nerve fiber branches include important landmarks for further disease study. It is, therefore, useful to detect all the branching points for quantitative analysis of nerve fibers. In this section, the corneal nerve fiber branches are located by employing the orientation response based branching detection in orientation scores [41]. An orientation score data representation is obtained by mapping a 2D image into a 3D space of positions and orientations R 2 S 1 via a wavelet-type transform. This transform is completed through convolution with an anisotropic wavelet ψ ∈ L 2 (R 2 ): where R θ = cos θ − sin θ sin θ cos θ denotes a counterclockwise rotation over angle θ . In this work, we choose cake wavelets [27] for ψ. Cake wavelets are quadrature filters. The real part shows the locally symmetric structures like curvilinear ridges/lines, while the imaginary part represents the antisymmetric structures like edges. For branching detection, we choose the real-valued wavelets to analyze the nerve fiber geometry. The orientation scores are functions defined on R 2 S 1 by lifting 2D images to the 3D coupled space of positions and orientations. When the image data is transformed to orientation scores, the anisotropic convolution kernels will perform the corresponding rotations and translations that follow the group operations defined in R 2 S 1 . Therefore, the rotating coordinate system and the rotating-invariant operations can always match each other for the processing of orientation scores.
In this work, we use 24 discrete orientations (from 0 to 2π) to map corneal nerve fibers into the orientation score representation. As shown in Fig. 4b, curvilinear structures like corneal nerve fibers are lifted into a 3D orientation space (x, θ) with a corresponding orientation response U f (x, θ) at each spatial location x. To avoid the false detection of topological outliers, the scores of dominant orientations are obtained by only considering the magnitude larger than a threshold value t = 1.2σ , where σ denotes the standard deviation of U f . In general, a corneal nerve fiber branching point splits into three separate segments at different directions. Thus, it can be detected by automatically counting the number of dominant orientations. In the orientation score domain, since we already have the responses U f (x, θ) of each pixel location at different orientations θ , we can efficiently obtain the number of dominant orientations by detecting the local maximum along the θ dimension, as shown in Fig. 4c, d. In principle, the proposed branching detection module also works on original images. However, the accuracy of the detection on original images can be affected by image noise and background artifacts. The low contrast of nerve fibers in original images results in relatively weak orientation responses, while the enhanced nerve fibers can produce high responses at branches. To better illustrate the improvement, in Fig. 5, we show the qualitative performance of branch detection on an original image and an enhanced image. We can observe that the orientation responses at branches of the original image are in general lower than the same branches in the enhanced nerve fiber images. The noise and background artifacts in the original image cause inaccurate responses due to unstable local maximums. Thus, the missing detections at low-contrast junctions can easily happen without enhancement.

Materials
We use the corneal nerve fiber images of Maastricht (CNM) dataset provided by The Maastricht Study [42], the Netherlands. This dataset includes 142 composite images (92 healthy individuals with normal glucose metabolism (NGM), 4 impaired fasting glucose, 16 impaired glucose tolerance, and 30 type-II diabetes (DM2) subjects) with a resolution of 1536 × 1536 pixels. In each of the images, the nerve fibers were manually annotated. The binary mask of each image was also manually created to define the region of interest.
We also validate both our segmentation and tortuosity measurement methods on the corneal nerve tortuosity (CNT) dataset [43] which consists of 30 images from the sub-basal corneal nerve plexus, acquired in normal (6 images, 20%) and pathological subjects, including diabetes (10 images, 33.3%), pseudoexfoliation syndrome (8 images, 26.7%), and keratoconus (6 images, 20%). This dataset was collected using the Heidelberg Retina Tomograph II with Rostock Corneal Module (HRTII32-RCM) confocal laser microscope (Heidelberg Engineering GmbH, Heidelberg, Germany). All the images consisting of 384 × 384 pixels cover an area of 400 × 400µm. The manual grading of the corneal nerve tortuosity is provided for this dataset.

Performance measurements
To measure the performance of corneal nerve fiber classification, we calculate the following metrics: true positives (TP), true negatives (TN), false negatives (FN) and false positives (FP). A pixel that is classified as nerve fiber in both the ground truth and the segmentation result belongs to TP, while a pixel that is classified as non-fiber in both the ground truth and the segmentation result is considered as TN. A pixel which is wrongly identified as non-fiber in the segmentation result belongs to FN, while a non-fiber pixel which is classified as a fiber pixel in the segmentation result is taken as FP.
The statistical performance measurements, including sensitivity (Se), specificity (Sp), and accuracy (Acc), are used to evaluate the global performance of a binary classification system. These measurements are given by where N = TN + TP + FN + FP. The Dice similarity coefficient (DSC) and the Matthews correlation coefficient (MCC) within the field of view (FOV) are reliable metrics for evaluating an unbalanced dataset with two classes of different sizes. In fact, this is the case in the corneal nerve fiber segmentation, where the percentage of nerve fiber pixels is only around 9-14%. The DSC only counts true positives in both the numerator and denominator, and the MCC is a correlation coefficient between the manual segmentation and the ground truth. The DSC and MCC are given by where S = (TP + FN)/N , P = (TP + FP)/N . The MCC value ranges from − 1 (completely incorrect classification) to 1 (perfect classification), and DSC from 0 (completely incorrect classification) to 1 (perfect classification). Moreover, the precision-recall curve (PRC) is used to evaluate the quality of the obtained nerve fiber probability maps since it is very suitable for measuring the performance of imbalanced datasets like corneal nerve fiber images. The PRC is formed by plotting the precision (TP/(TP + FP)) versus the recall (TP/(TP + FN)) with respect to the varying threshold value t. All the performance measurements are taken within the FOV.

Corneal never fiber segmentation on the CNM images and challenging structures
In Fig. 6, we show qualitative corneal nerve fiber segmentation results on the CNM dataset using our proposed method. We can see that our method achieves very good segmentation performance on extracting corneal nerve fibers. Our method is able to preserve tiny structures with low contrast and high illumination changes, as shown in the small patches in column 3-4 of Fig. 6. Even small disconnections at low-visible regions are recovered after enhancement. Although the corneal nerve fiber images exhibit background artifacts and high-noise characteristics, we can still observe from Fig. 6 that those factors have limited influences to the segmentation performance. This is due to that most of the false positives actually come from the missed annotation of weak nerve fiber structures in the ground truth images. This indicates that the proposed FCN framework is capable of preserving the corneal nerve fiber structures and meanwhile suppressing the background noise. In Fig. 7, several examples of complex nerve fiber structures are given to show the enhancement ability of the proposed method. In rows 1 and 3 of Fig. 7, we can observe that our method is better at dealing structures with strong occlusions. The background noise and non-fiber artifacts are reduced after enhancement. In row 2, we see that tiny nerve fibers with disconnections are better highlighted when compared to the results from the other two methods.  [44] and LIDOS [29], respectively This shows that the proposed deep FCN-based segmentation method produces high-quality nerve fiber images for subsequent analysis.

Quantitative evaluation of the segmentation performance
To investigate the influence of our proposed corneal nerve fiber segmentation method, we compare our segmentation results on the 142 corneal nerve fiber images with the results obtained from the state-of-the-art segmentation approaches [29,44] as shown in Table 1. The sensitivity, specificity, DSC and MCC metrics are used to evaluate the global performance since they are suitable for evaluating the corneal nerve fiber images which have imbalanced classes. For the sake of equal comparison, we obtain the best performance of the Frangi vesselness and LIDOS (left-invariant derivative filters in orientation scores) filters with the same preprocessing steps as in our approach. We can see that the proposed method achieves significantly higher DSC and MCC values than the other methods. All the three methods obtain very high specificities due to the large portion of background pixels in the corneal nerve fiber images. Nevertheless, the proposed method still achieves the highest sensitivity of 0.8632, DSC of 0.8433 and MCC of 0.8431 at the highest specificity of 0.9978 compared to the other two methods.
The receiver operating characteristics (ROC) curve is often used to evaluate the quality of a binary segmentation. However, in this work we employ the precision-recall curve (PRC) as an alternative since it is more suitable to represent the segmentation performance of the extremely imbalanced corneal nerve fiber images. Figure 8 shows the precision-recall curve of the corneal nerve fiber segmentation results. We can see that the proposed method gives much better and stable performance in the low threshold value region (right part of Fig. 8).

Fig. 9
Examples of our corneal nerve fiber segmentation and tortuosity measurement on the CNT dataset To better validate the proposed nerve fiber segmentation and tortuosity measurement techniques, the CNT dataset is processed to evaluate the performance of our pipeline. Since this dataset only provides the tortuosity grading of each image without having any manual annotation of the nerve fibers available, we employ a crossing-training to obtain the segmentation results using the classifier trained from the CNM dataset. Figure 9 shows the qualitative performance of our segmentation method on this dataset. We can observe that most of the fibers are perfectly extracted from those CNT images based on our trained classifier from the CNM dataset. Afterwards, we calculate the tortuosity measure κ μ of each image in the CNT dataset using our tortuosity measurement method described in Sect. 2.3. In Fig. 9, we can see examples of showing the measured fiber tortuosity values 0.0938, 0.2108 and 0.8539, respectively, on three typical images from different stages labeled by experts. The tortuosity measurements intuitively show the fiber tortuosity variations of different cases, which indicate the reliability of our tortuosity calculation method.

Validation of the branching biomarker detection
The orientation score-based approach described in Sect. 2.4 provides us with an efficient way for localizing corneal nerve fiber branches. Based on the FCN-enhanced nerve fiber map, the orientation responses at each branching location are better optimized under the well-preserved curvilinear continuity. Unlike other curvilinear structures in medical images, the corneal nerve fiber junctions in general only include branching structures (with three directions) but not crossings (with four directions). Thus, in our setting we mainly consider the detection of critical location with three maximum responses. We evaluated the nerve fiber branching detection approach on both the CNM dataset and the CNT dataset. To quantitatively evaluate the detection performance, we manually labeled all branching points of 20 images from each dataset for comparing with the automatic detection approach. The proposed pipeline successfully achieves high accuracy of 98.3% and 99.5% on the CNM and CNT images, respectively. In Fig. 10, we provide typical examples from both datasets to qualitatively present the performance of branching detection on corneal nerve fibers. Figure 10a shows the result on an image from the CNM dataset and Fig. 10b-e shows the results on four images from the CNT dataset. We can observe that almost all the branching points in these images are accurately located using the orientation score responses based method.

Conclusion
In this work, we proposed a fully automated solution for corneal nerve fiber touristy measurement and branching detection in microscope images. In the first phase, an automatic deep learning-based framework was proposed to enhance and extract corneal nerve fiber images. We applied intensity normalization to effectively address the homogeneity of the images. A FCN was proposed to address the challenges of low nerve fiber visibility and high-background artifacts. The proposed framework achieves promising quantitative and qualitative results based on the evaluation on two datasets. In the second phase, an automatic tortuosity calculation and an accurate branching detection methods were established to extract valuable geometric and topological biomarkers from corneal nerve fiber images. Both methods have been validated on the CNM and CNT datasets to show their effectiveness and reliability. In the future, we will use those extracted biomarkers to analyze clinical data for disease discrimination.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.