Meibomian glands (MG) are special sebaceous gland located on the tarsal plate of eye (Erich et al., 2011). They produce meibum to prevent tear from over evaporation (Nichols et al., 2011), maintain the surface tension of the tear film, and trap tears between its oily edge and the eyeball. Healthy MG are elongated in shape, arrange in parallel and throughout the entire tarsal plate (Djalilian, 2018). Functional and/or structural problem of MG may cause meibomian gland dysfunction (MGD) (Driver and Lemp, 1996). MGD is a chronic, diffuse abnormality of the meibomian glands, commonly characterized by terminal duct obstruction and/or qualitative/quantitative changes in the glandular secretion (Den et al., 2011; Nichols et al., 2011). MGD often leads to tear film alteration, ocular surface disease, and is the leading cause of evaporative dry eye (Turgut et al., 2018), which can seriously affect the patient’s life and increase the global public health and financial burden (McDonald et al., 2016).

Meibomian gland area loss is an important index to evaluate MGD (Arita et al., 2014). MGD can be directly observed using meibography, which is an optical imaging technique allowing visualizing MG morphology in vivo (Pult and Nichols, 2012). Recently, many non-contact infrared meibography methods were developed, making the MGD diagnosis process more patient-friendly and less time-consuming (Pult and Riede-Pult, 2012; Wong et al., 2019; Xiao et al., 2020; Hwang et al., 2021). These devices allow users to capture high-resolution images of meibomian glands in a short period of time, which can provide sufficient experimental material for MG analysis. Quantification of the area of meibomian glands loss is of importance when assessing MGD. To date, automatic methods based on image processing techniques have been developed for the automated assessment and classification of MGD (Arita et al., 2014; Celik et al., 2013; Koprowski et al., 2016; Koprowski et al., 2017; Liang et al., 2017; Llorens-Quintana et al., 2019; Pult and Riede-Pult, 2013; Xiao et al., 2021). In recent years, MGD automatic analysis methods based on convolutional neural network (CNN) have been developed rapidly. These works automatically segment the eyelid and MG, and calculate the loss rate and analyze the morphological parameters of MG (Wang et al., 2019). In order to promote the segmentation performance of MG, contrast limited adaptive histogram equalization (CLAHE) will be used as an image preprocessing step to enhance the contrast of MG images (Prabhu et al., 2020; Dai et al., 2021).

However, the MG images used in these works were acquired from a single device. The preset parameters of the method or the trained model are for the specific data domain used in the experiment. When testing these methods with images from other distribution domains (such as different image modalities, or different acquisition devices), a lower performance is usually obtained (Yan et al., 2019). Such problem is called “distribution shift” (Jo and Bengio, 2017). To reduce the performance gap, an effective method is to reduce the distribution domain distance of the data. Training a generative adversarial network (GAN) to generate fake images between two different domains is a common method to improve the generalization of the model to data from different domains (Perone et al., 2019; Yan et al., 2019); but it requires a large number of training samples. In some medical image processing tasks, histogram specification (HS) has been reported as an effective method for contrast enhancement and improving model performance on images of different modalities, such as between MRI and CT images (Naseem et al., 2021).

In this study, we focus on two open questions in the field of automated meibomian gland analysis:

1) In the CNN-based MG analysis method, does CLAHE improve the network performance?

This study aims at developing and validating an MG segmentation method based on CNN. Whether HS and CLAHE as preprocessing methods can effectively improve the segmentation performance of the CNN model was also investigated in this study. MG images captured from different devices were divided into different testing sets based on different preprocessing methods, and then segmented by CNN model to calculate MG loss rate. All the results will be compared with the ground-truth by the clinicians.

Histogram specification (HS) enhances the brightness and contrast of the input image and transforms the input image into an image with a similar shape to the template histogram. HS has the advantage of simplicity and low computational cost (Xiao et al., 2018). HS is very common in the preprocessing stage of medical images. HS is used to match the histogram of the input image to the template image to initialize and avoid gradient explosion before using CNN to segment the kidney based on CT images (da Cruz et al., 2020). Naseem et al. (2021) based on HS to enhance images of different modalities. They enhance low-contrast CT images based on MRI images based on the second-order distribution.

The raw MG image often has low contrast. In other word, the pixel intensity of glands and background have little difference, as a result, it is difficult for the observer to separate the glands from background, especially when the image is blur. CLAHE is applied for the meibography image enhancement. CLAHE enhances an image by limiting the height of the local region histogram, such as a 64 pixels neighborhood block, thereby overcome the global uneven illumination and noise amplification problem. CLAHE is often used as a preprocessing method to enhance the performance of MG segmentation (Prabhu et al., 2020; Dai et al., 2021).

This study involves 287 subjects (age: 56.1 ± 17.2 years old, 83 men and 204 women) collected by Beijing Tongren Hospital. The purpose and possible consequences of the study were explained to all involved subjects. Exclusion criteria included 1) ocular allergies, 2) history of ocular surgery, 3) history of ocular trauma, 4) other eye diseases, 5) long term or frequent contact lens use and 6) Any other eye or systemic disease known to affect the tear film. Excluded images include: 1) Images included something other than the eyelids and their surrounding tissueswere, 2) Images were not sufficiently clear for automatic analysis were excluded. 3) Patients’ eyes exhibited excessive meibomian lipid secretion. The study was approved by the Ethical Committee of the Beijing Tongren Hospital and was conducted in accordance with the tenets of the Declaration of Helsinki. A total of 1,074 images were collected by professional clinicians, including 888 MG images collected from Keratograph 5M (K5M; OCULUS Optikgeräte GmbH, Wetzlar, Germany) and 186 from KANGHUA DED-1L (KH). Note that I k 5 m is used to represent the image set collected from K5M, and I k h is used to represent the image collected from KH. Then, I k 5 m was divided into two sets, 648 for training set used for neural networks training, and 240 for testing set to appraise the algorithm performance. And I k h is used as an external testing set to test the performance of our models on data from different device. It should be noted that as a retrospective study, some of the 287 patients have received treatment and collected follow-up data at different times, therefore there are more meibomian gland images than the subject number.

The ground-truth annotation was completed by a senior clinician using imageLabeler application in MATLAB (9.5.0.944444, R2018b, Java 1.8.0_152-b16). For labeling the eyelid region, the upper edge is defined at the opening of the gland, the lower edge is defined at the edge of proximal tarsal plate, and the horizontal borders is defined at the top and bottom borders intersected (Wang et al., 2019). Value 1 and 2 are used to denote the eyelid region and the MG region, respectively. Examples are shown in Figure 1.

Two models were trained with the training set, including: 1) the model was trained with the k5m generation graph; 2) K5M images enhanced by CLAHE were used for training. Each training image and its corresponding mask were resized to 256 × 256 and augmented by gamma transformation, rotation, blur, noise addition, and image flip. The data augmentation increased the training images by 12 times. Training and evaluating of the proposed models were performed on a computer with an Nvidia GeForce GTX 3090 GPU. The deep learning models were implemented based on PyTorch (version 1.1.0) package in Python. All models were trained using the Adam optimizer (α = 0.9, β = 0.999), with an initial learning rate of 0.0003 and decays 0.8 times every five epochs. The batch size is set to 32 and the maximum epoch number of 50. L1 regulation was applied to prevent over-fitting (Hawkins, 2004). The loss function used in the experiment is Dice loss function (Jadon, 2020), which is represented by Eq. 2. L D i c e = 1 − 2 · | y ∩ y ′ | | y | + | y ′ | (2)

To explore the influence of HS and CLAHE on model performance, k5m images and KH images were divided into different testing sets, which were: 1) I k 5 m and I k h ; 2) apply HS to KH image with I k 5 m as the guide images; 3) preprocessed I k 5 m and I k h using CLAHE. The specific division is shown in Table 1.

The performance of the MG segmentation models was assessed by dice similarity coefficient (DSC), recall, and precision. These results were reported in Table 1 and Table 2. The results in Table 1 show that all images have close results. The maximum value of each indicator is marked in bold. Table 3 shows the comparison between segmentation results of different preprocessing methods on two image set. Note that there is no statistical difference in all different testing sets, except that the precision between the I k h − C L A H E and I k h − H S − C L A H E is statistically different. Similar results can also be found in the I k 5 m and I k 5 m − C L A H E , with no statistically significant difference between the evaluation results of the segmentation model regardless of whether the images were processed with CLAHE or not.

Figure 4 shows some visualization of MG segmentation results. Most of the images in both the internal testing set and the external testing set have segmentation results with high DSC score, and the network predicted segmentation masks and ground-truth masks have high similarity and coincidence. But some images in the external testing set have very low segmentation results, as shown in Figure 5, both the eyelid region and the MG region are under-segmented, with rough edges and anomalous masks. These under-segmented images have dissimilar features when compared with the predicted results on average level.

Figure 6 and Figure 7 show the direct comparison results between the MG loss rate and GT for the external and internal testing sets. Most of the predicted MG loss rates are distributed near the ideal line (predicted MG loss rate equals to ground-truth MG loss rate), which indicates that the MG loss rates calculated according to the predicted segmentation results of the network are relatively close to the ground-truth. Root Mean Squared Error (RMSE) in Table 4 was used to quantify and compare the MG loss rate with the ground-truth. Testing set I k 5 m − C L A H E has the smallest RMSE when compared with ground-truth (RMSE = 0.09).

Meibomian gland morphology analysis is an important index to evaluate MGD. However, the current manual grading methods have the problems of heavy workload, low efficiency and large error, which is not conducive to the standardized diagnosis and treatment of MGD and dry eye. In this study, we developed and evaluated the automatic segmentation method of the eyelid area and the MG area based on CNN and automatically calculated MG loss rate. This method is evaluated in the internal and external testing sets from the two meibography devices. In addition, we also tested that the pre-processing method of HS and CLAHE as MG images has not improved significantly.

In this experiment, the segmentation results of the external testing set ( I k h , I k h − H S , I k h − C L A H E , and I k h − H S − C L A H E ) from another device are lower than the internal testing set ( I k 5 m and I k 5 m − C L A H E ). Such performance gap is common in segmentation tasks based on CNN. Even though data augmentation was used in the training phase, data augmentation in the distribution domain of the training set ( I k 5 m and I k 5 m − C L A H E ) cannot close the distribution gap, because domain shift is systematic. As shown in the histogram in Figure 8, I k h and I k 5 m have very different pixel distributions. These image differences lead to a performance gap in the segmentation performance of the model between the two different testing sets. We used HS to make the histograms much more similar, but the obtained network segmentation results did not improve significantly. t-SNE (Van der Maaten and Hinton, 2008) was used to gain insight into the difference between the two sets of images before and after using HS, as shown in Figure 9A. The images are applied dimensionality reduction and projected onto a 2D plane. There is a large gap between the distributions of I k 5 m and I k h . Also, HS does not effectively reduce the distribution gap, although the histograms had been matched. Changing the low-dimensional image feature of gray distribution will not affect CNN to extract and learn the high-dimensional semantic information of the image.

CLAHE can effectively enhance the visual effect of MG images, which is beneficial for doctors to read and measure MG images manually. However, CLAHE does not significantly improve the segmentation results of neural networks, because CLAHE also fails to close the distribution gap between the distributions of I k 5 m and I k h , as shown in Figure 9B. CNN have been shown to be feasible in medical image segmentation tasks, but when used in practice, network models can suffer huge performance degradations if image data from other distributed domains are involved. Such performance gap can be explained using the independent and identically distributed (i.i.d.) assumption of statistical learning: a network model that is well-trained in the source distribution domain does not necessarily achieve similar high performance on datasets with the same distribution as the source distribution domain performance (Zhu et al., 2017). Adjusting the pixel distribution of images by HS and CLAHE did not improve the prediction results of the network model, although the images were visually very similar to the source distribution domain. Through tsne-2D visualization, it can be found that HS and CLAHE do not close the distribution gap.

In this study, we developed and validated a CNN-based automatic MG analysis method based on MG images acquired by two different devices, K5M and KH. Predictions of MG loss rates were in high agreement with the gold standard obtained by physicians. We also found that in this study HS and CLAHE do not significantly improve the performance of CNN for segmenting MG, which may indicate that in future deep learning-based MG analysis tasks, no additional computational cost is required for image preprocessing phase.