Image-based quantitative analysis of tear film lipid layer thickness for meibomian gland evaluation

Dry eye syndrome is one of the most common ocular diseases, and meibomian gland dysfunction (MGD) is the leading cause of evaporative dry eye syndrome. When the tear film lipid layer becomes thin due to obstructive or hyposecretory meibomian gland dysfunction, the excessive evaporation of the aqueous layer can occur, and this causes evaporative dry eye syndrome. Thus, measuring the lipid layer thickness (LLT) is essential for accurate diagnosis and proper treatment of evaporative dry eye syndrome. We used a white LED panel with a slit lamp microscope to obtain videos of the lipid layer interference patterns on the cornea. To quantitatively analyze the LLT from interference colors, we developed a novel algorithm that can automatically perform the following processes on an image frame: determining the radius of the iris, locating the center of the pupil, defining region of interest (ROI), tracking the ROI, compensating for the color of iris and illumination, and producing comprehensive analysis output. A group of dry eye syndrome patients with hyposecretory MGD, dry eye syndrome without MGD, hypersecretory MGD, and healthy volunteers were recruited. Their LLTs were analyzed and statistical information—mean and standard deviation, the relative frequency of LLT at each time point, and graphical LLT visualization—were produced. Using our algorithm, we processed the lipid layer interference pattern and automatically analyzed the LLT distribution of images from patients. The LLT of hyposecretory MGD was thinner (45.2 ± 11.6 nm) than that of dry eye syndrome without MGD (69.0 ± 9.4 nm) and healthy volunteers (68.3 ± 13.7 nm) while the LLT of hypersecretory MGD was thicker (93.5 ± 12.6 nm) than that of dry eye syndrome without MGD. Patients’ LLTs were statistically analyzed over time, visualized with 3D surface plots, and displayed using 3D scatter plots of image pixel data for comprehensive assessment. We developed an image-based algorithm for quantitative measurement as well as statistical analysis of the LLT despite fluctuation and eye movement. This pilot study demonstrates that the quantitative LLT analysis of patients is consistent with the functions of meibomian glands clinically evaluated by an ophthalmologist. This approach is a significant step forward in developing a fully automated instrument for evaluating dry eye syndrome and for providing proper guidance of treatment.


Background
Dry eye syndrome (DES) is one of the most common ocular diseases affecting more than 5% of the US adult population [1]. Dry eye syndrome is composed of aqueous deficiency DES and evaporative-type DES. The most common cause of the evaporative-type DES is the meibomian gland dysfunction (MGD) [1][2][3]. In the eye, the preocular tear film comprises of mucous, aqueous, and lipid layers. The lipid layer, formed by oil (meibum) secreted from meibomian glands, prevents excessive evaporation of the aqueous layer. In obstructive or hyposecretory meibomian gland dysfunction, thinning of the lipid layer leads to excessive evaporation and develops evaporative-type DES [4,5]. Recently, numerous studies have been conducted on the correlation between DES and the lipid layer thickness (LLT) [6][7][8][9][10]. Therefore, accurate measurement of LLT has become more important.
In normal condition, the LLT has been estimated and reported to be approximately 70 nm in the eye [11,12]. The lipid layer thickness can be estimated by analyzing the color intensity patterns generated by the optical interference from multiple reflections at the interfaces of air-lipid and lipid-aqueous layers [13]. Various methods utilizing interference patterns have been used to characterize the lipid layer thickness: gooseneck light [14], bio differential microscope [15], Polaroid filters [16], monochromatic light [17], spectral-discrimination [18], and a simple interferometer made of paper tool for lipid layer evaluation [19]. Several tear interference imaging devices have also been developed, such as the DR-1 tear interference camera (Kowa Co., Nagoya, Japan) [20,21], LipiView interferometer (TearScience Inc, Morrisville, NC) [9,22] and Lipiscanner 1.0 (Visual Optics, Chuncheon, Korea), a cost-effective add-on to an existing slit lamp biomicroscope. So far, only the LipiView device can provide quantitative values of the lipid layer thickness.
While the DR-1 can qualitatively visualize the interference pattern of lipid layer [23], the LipiView interferometer can quantitatively measure the average lipid layer thickness. The main proprietary algorithms and high-speed computers in these systems capture the reflected color from lipid layer at a rate of approximately 14 million pixels per second to complete the evaluation. The spatially modulated light source allows the elimination of unnecessary background images and stray light. The processed output is expressed as interference color that correlates with the thickness of the lipid layer [24].
However, complicated interferometry systems are costly (Lipiview) or may provide qualitative analysis (DR-1). By utilizing a recently developed low-cost custom-made Lipiscanner 1.0 system for quantitative measurements, we developed an image-based algorithm to automatically define the region of interest (ROI) on the iris and to evaluate the thickness of the lipid layer, even with irregular pupil movements during data acquisition. (see Additional file 1: Movie, Additional file 2: Movie, Additional file 3: Movie, Additional file 4: Movie.) Also, we performed color compensation in our algorithm to obtain a precise measurement of the lipid layer thickness that is based on the Fresnel equation [25]. To validate our algorithm, we performed a feasibility analysis on the video data obtained from twenty patients and fourteen healthy volunteer who were affected by either hyposecretory MGD, dry eye syndrome without MGD, or hypersecretory MGD.

Theory
A colorful interference pattern is caused when the light is reflected from the top and bottom boundaries of a thin film such as an oil film floating on a bubble or water. Interferences by thin films display different colors depending on the thickness of the film, from dark color caused by thinner film area to brighter color caused by thicker film area. Using Snell's law, where OPD means "optical path difference, " n m means "refractive index of meibum, " d means "lipid layer thicknesses" and β means "refraction angle. " The phase difference ( ) between the two interfering lights can be calculated using the wavelength of light ( ) and OPD as follows: Meanwhile, the intensity of light is calculated by the Fresnel equations. The Fresnel equations are a pair of equations that describe the reflectance and transmittance of a surface between two media having different refractive indices. The reflectance and the transmittance are coefficient values representing the ratio of incident light that is reflected or transmitted. According to the Fresnel equations, the ratio of the reflected and transmitted wave's complex electric field amplitude to that of the incident wave (reflection and transmission coefficients) for s-and p-polarization are r s , r p , t s and t p : where θ 1 and θ 2 mean incident angle and refractive angle, n 1 and n 2 indicate refractive index of air and meibum.
Then, the reflectance and transmittance are as follows: where R and T represent the reflectance and the transmittance, respectively.
(3) r s = n 1 cos θ 1 − n 2 cos θ 2 n 1 cos θ 1 + n 2 cos θ 2 (4) t s = 2n 1 cos θ 1 n 1 cos θ 1 + n 2 cos θ 2 (5) r p = n 2 cos θ 1 − n 1 cos θ 2 n 2 cos θ 1 + n 1 cos θ 2 (6) t p = 2n 1 cos θ 1 n 1 cos θ 2 + n 2 cos θ 1 With these coefficients, all intensities of rays can be defined in a simple reflectiontransmission model. In this model, we defined R m (m, a nonnegative integer) to be the ray's wave function which is reflected m times between the lipid-tear interface and the air-lipid interface. R 0 is the wave function of the ray which is reflected directly at the surface of the lipid layer. Hence, where I R n (nonnegative integer n) is the intensity of the light R n , and R i is the incident light's wave function. R ij is the coefficient of the reflectance from i layer to j layer (1-air, 2-lipid layer, 3-tear layer) and T ij corresponds the coefficient of the transmittance.
In this model, OPD is present at every new reflection-transmission process. Therefore, the phases of the rays are described as follows. Now we have each ray's phase and intensity. From these, the interference ray can be calculated in the complex number plane as follows.
I INT ( , d) can be obtained from the amplitude of the wave function R INT ( , d), which is used to obtain the RGB value as follows: Here, I INT ( , d) is the intensity of interference light obtained with the wavelength ( ) and lipid layer thickness (d), and the functions R STDOBS , G STDOBS , B STDOBS of represents the value of the color whose wavelength corresponds to in the CIE1964 RGB standard observer. The wavelength range is from 360 to 830 nm in the visible light range provided by the CIE 1964 standard observer.
In this way, we can calculate the interference color as RGB values at different thicknesses of the meibum in the range of 0-240 nm. (9)

Methods
We developed our algorithm to measure the lipid layer thickness quantitatively by analyzing the lipid layer interference patterns on tear film with the Lipiscanner 1.0 which consists of a LED panel (115 mm × 58 mm) covered with a polycarbonate diffuser for homogeneous illumination (Fig. 1). This simple device was used to observe the lipid layer of the tear film with a slit lamp microscope for ophthalmology and a scientific complementary metal-oxide semiconductor (sCMOS) camera. A patient's head is placed in a fixed position on the head-chin rest of the slit lamp microscope and white light from the LED is irradiated onto the lipid layer of the eye. The white LEDs provide a color temperature of ~ 6500 K. The color from the white light interference is used to assigning the LLTs to ROI image pixels. The measurements presented in this paper represent the LLT of the precorneal tear film in the inferior iris region which was illuminated by white light.

Instrument setup
We used the Lipiscanner 1.0 to visualize the lipid layer by white light interference, and the captured videos were used for measuring the lipid layer thickness (Fig. 1). In this setup, the LED light was illuminated onto the inferior cornea of the patients, which was then reflected onto the slit lamp microscope. The illuminating beam was not focused to a point but was spread over parts of the cornea by a diffuser. An sCMOS camera (Guppy Pro F-503, Allied Vision) was used to record the video of the interference pattern of lipid layer within the cornea region.

Patients with dry eye syndrome
This study followed the tenets of the Declaration of Helsinki and was approved by the Institutional Review Board of Chuncheon Sacred Heart Hospital. Twenty patients with dry eyes and fourteen healthy volunteers were included in this pilot study: six patients with hyposecretory MGD (low delivery of meibum) (patient group I), seven patients with dry eye syndrome without MGD (patient group II), seven patients with hypersecretory MGD (high delivery of meibum) (patient group III), and fourteen healthy volunteers as a control group (patient group IV). An ophthalmologist who is also a meibomian gland expert (Dr. Ho Sik Hwang) categorized the patients into groups I, II and III after Fig. 1 Schematic diagram of the optical system for thin film interference measurement system evaluating their medical histories [26], slit lamp examinations to measure tear break-up time [27], corneal stain, lid margin abnormality, meibum volume and quantity, Schirmer's test [28] for tear production measurement, and meibography for morphological evaluation of the meibomian gland.

Image processing
One of the main challenges in the processing of the eye images in our setup was that the pupil often moves suddenly during the image acquisition. To address this, we developed an algorithm that was robust even when the pupil location changed. For each patient, our algorithm processed an image sequence of 2.5 s duration (a subset of the original sequence, i.e. 75 images), and we discarded images where the pupil was not visible or was partially occluded (e.g. caused by blinking). Afterward, we started with the darkest spot in the image (which is a point in the pupil) and performed a region-growing process to extract the pupil region. We repeated this process with a different starting point to extract the iris region, from which we could extract our region of interest (ROI) of where the interferometry colors existed. Finally, we compensated colors that could be altered by the illumination light or the ambient room light so that accurate measurements of the lipid layer thickness could be achieved. The entire image processing part of the algorithm to measure the thickness of the lipid layer runs through six phases, as shown in Fig. 2.
The Additional file 1: Movie S1 shows the procedure of image processing algorithm.

Phase 1: Exclude unnecessary frames
First, we have to select the image frames in the raw video for analysis and discard frames that are not suitable for use. The resulting extracted video must satisfy the following conditions: The center of the pupil in each frame must be near the center of the screen (to increase recognition accuracy), and there should be no change in image brightness and zoom level between image frames. We achieve this by filtering out frames that have no ROI, including cases when the eye closes. When eye blinking occurs, the light from the LED is strongly reflected from the eyelid, and the entire image becomes brighter. We define B opened as the brightness when the ROI is visible, and B closed as the brightness when the eye closes. Then the relationship of the values is B opened < B closed , where B opened is set to the mode of full frame brightness data, and B closed is set to the highest value of all brightness data. The threshold value to filter out frames including eye blinking is B opened + 0.33 B closed − B opened . The coefficient of 0.33 helps to filter out the unnecessary frames and prevents false positive frames from being included in the ROI data (see Additional file 5: Image S1).

Phase 2: Find the pupil and iris regions
Under our current image acquisition conditions, the darkest region (21 × 21 pixels) of an image in our sequence is located inside the pupil that can be used to demarcate the pupil region. Starting with the centroid of the darkest region in the image, we apply the flood-fill algorithm with 8-directions to select all pixels that belong to the pupil [29]. The Flood-fill algorithm is an algorithm that starts at a point and selects a connected group of pixels where the color (or brightness) difference between the examined pixel and the starting pixel is smaller than some selected threshold value (see Additional file 6: Image S2). Applying the flood-fill algorithm for extracting the pupil, the boundaries of the region may not be smooth and contain a portion of the iris. Thus, it is necessary to remove unwanted iris region as it affects the centroid by blurring the image followed by re-applying flood-fill algorithm with different threshold value (see Additional file 7: Image S3).
Once the image has been blurred, additional flood-fill algorithm with a lower threshold value can be used to capture the shape of the pupil with smooth boundaries without including the iris. After the pupil region is successfully extracted, the centroid of this region can be calculated which gives an approximation to the location of the center of the eye. If the location of the center of the eye obtained after correction using blurring image and the flood-fill algorithm is far from the darkest region in the pupil, we take the centroid of the darkest region as the center of the eye instead while we may partially lose the ROI.
To extract the iris region, we apply the same flood-fill algorithm again with a starting point outside the pupil (i.e. a point in the iris). We then perform the Canny edge detection [30] on the resulting region to extract the boundary pixels. However, as the iris is partially covered by the eyelid most of the time, the resulting boundary consists of many false boundaries of the iris. To overcome this, we examine only the boundary pixels that have almost vertical edge orientations and are within an empirically pre-defined distance from the detected center of the eye (in our data, the diameter of iris was approximately 480 pixels. However, this number depends on the camera zoom state and patient's eye trait. Thus, in practice this iris size has to be determined using the iris from image.) It is because in order to get as much ROI as possible, we need the exact radius of the iris but already know the approximate value. From these boundary pixels, we can take their average distances from the center of the eye and estimate the radius of the iris.

Phase 3: Extract the region of interest
The ROI should be the region within the iris that shows interference colors that should appear in the bottom half of the iris image. Thus, we choose our ROI to the area within 80% of the measured iris to prevent the lower eyelid from being included in the ROI and to exclude the sclera as part of ROI in case the centroid of pupil is inaccurate. We then crop away the regions above the center of the circle, the regions outside of the circle, and for the region within the circle, we also crop away small regions on the left and right that are farther than 80% of the radius of the new inner circle. This is to prevent the appearance of saturated pixels from the white sclera of the eye. After that, pixels that have a brightness less than the average brightness of the region are also removed. The resulting region is our ROI. Figure 3 shows some of the intermediate images in phases 2 and 3 and the final ROI.

Phase 4: Subtract iris color from region of interest
The colors in the ROI originate from a combination of the white light interference and the iris color. We need to subtract the iris color component from the ROI, or otherwise, it can affect the subsequent color analysis. We achieve this by finding the color of the iris in the phase 3 circle outside of the ROI, and then subtract it from the ROI.

Phase 5: Correct for illumination colors
Depending on the color and brightness of the room lighting or camera exposure value at the time of image acquisition, the colors in the image may appear to be different from the actual ones. Also, the color temperature of the white light that we used to calculate the theoretical color corresponding to the thickness is different from that of the LED light of the Lipiscanner. To compensate for this, we estimated how much the RGB values were biased from the white color of the sclera of the patient's eye captured under the same conditions and applied the corresponding correction to the ROI.

Phase 6: Assign thickness to all ROI pixels
We map each pixel's color to a thickness value of the lipid layer using a lookup table represented by the three-dimensional solid curve in the RGB space (Fig. 6). The lookup table is obtained by applying the Fresnel equations to the reflection-transmission model and then plotting the color output against a different lipid film thickness input. Color values that do not match any values in the lookup table are approximated by the color in the lookup table that has the closest Euclidean distance from it. In this way, the distribution of LLT variation within each frame and across image frames can be assessed, and the mean and standard deviation of LLT can be calculated.

Region of interest tracking
The automated ROI tracking technique and some intermediate results from Phase 2 and 3 of our algorithm are shown in Fig. 3 with one exemplary image frame.
• Preparing original images from a video clip (Fig. 3a).
• Searching for the darkest point in the image (the white dot inside the pupil) (Fig. 3b).
• Identifying of the pupil region of the eye (white color) (Fig. 3c).
• Determining the pupil of the eye (marked with red boundary) (Fig. 3e).
• Locating the centroid of the pupil (the red dot) (Fig. 3f ).
• Finding the boundary of the iris region (Fig. 3i).
• Determining the radius of the iris and drawing a circle (Fig. 3j).
• Defining a smaller subset region of the iris (Fig. 3k).
• Setting the radius of a smaller circle (Fig. 3l).
• Cropping the area above the center (Fig. 3m).
• Selecting the area inside of the small circle (Fig. 3n).
Setting the region of interest (Fig. 3o).

Statistical analysis of lipid layer thickness over time
The quantitative measurements of the lipid layer thickness from three representative patients in each group are summarized in Fig. 4. Patient I has hyposecretory MGD, patient II has dry eye without MGD, and patient III has hypersecretory MGD. The relative frequency distributions of the lipid layer thickness at selected time points in the three patients are shown in Fig. 4a-c. The average lipid layer thickness and the standard deviation per frame are also computed and plotted for the entire image sequence. The mean lipid layer thicknesses for patient I, II and III were measured as 35.6 ± 14.3, 75.2 ± 27.8, 120.3 ± 54.2 nm (average ± standard deviation), respectively from the data in Fig. 4d-f. A real-time display of the analysis of the lipid layer thickness for the three patients including means, standard deviations, and a relative frequency graph for each frame are included in the Additional files (see Additional file 2: Movies S2, Additional file 3: Movie S3, Additional file 4: Movie S4).

Visualization of lipid layer thickness information in space
We constructed a three-dimensional (3D) visualization of the structure of each patient's lipid layer for individual assessment of the lipid layer thickness of the patients' data shown in Fig. 5. This visualization with color coding for the lipid layer thickness can help us to examine overall LLT distribution and parts that are exceptionally thin or thick. The closer the ROI is to the lower eyelid, the thicker the LLT (the direction in which y increases). Due to the eyelash shadows on the ROI, we have a pixel that does not show a uniform LLT on the 3D surface plots. The LLT value in eyelash shadows is close to zero and is not included in the analysis. Figure 6 shows plots of the compensated image pixel's color data from the three patients. The color scale curve-3D representation of color bar on the right side-represents theoretical LLT values in three-dimensional RGB space. The individual dots of scatter plot represent the image pixel's RGB color data, and the color of each data represents the assignment of LLT based on the closest Euclidean distance from the color lookup table. Patient I image pixel data is distributed around the low LLT, Patient II data is near 100 nm, and Patient III data is distributed near the high LLT graph. The number of image data pixels in a graph is about 3000-5000. Data were obtained from three image frames, including the front and back frames of the image frame used to draw the 3D surface plots (Fig. 5). We downsized the number of data in the graph to one-third and used a scatter plot.

Comparison of the LLT and clinical data of patients with dry eye syndrome
As a proof-of-principle, we used a total of thirty-four patients with dry eye syndrome for this feasibility study. Their data were blindly analyzed with our image processing algorithm without any clinical information of the patients. Based on routine clinical examinations, an ophthalmologist categorized the patients into four groups: hyposecretory MGD (low delivery of meibum, the patient group I), dry eye syndrome without MGD (patient group II), hypersecretory MGD (high delivery of meibum, patient group III), and Group IV (healthy volunteers as a control group). The results of the LLT measurements for each category is compared in Fig. 7. Because changes in the LLT as captured in the patient images are slow and continuous, using a camera system at 30 fps was sufficient to obtain the results. With the limited data set, the LLT of hyposecretory MGD was thinner (45.2 ± 11.6 nm) than that of dry eye syndrome without MGD (69.0 ± 9.43 nm), and the LLT of hypersecretory MGD was thicker (93.5 ± 12.4 nm) than that of dry eye syndrome without MGD and healthy volunteers (68.3 ± 13.7 nm). There were statistical differences between the three groups (a one-way ANOVA).

Discussion
The purpose of this study was to develop an image-based algorithm for quantitative analysis of tear film lipid layer thickness (see Additional file 8: Image S8). Our algorithm was robust enough to work despite fluctuations and eye movements during the video recording. For example, the ambient room light did not affect our measurement as the illumination was relatively bright. In addition, we used a low-grade scientific camera with relatively high noise that has not been problematic.
We created a video clip to explain our algorithm involving the steps and results (see Additional file 1: Movie S1). We analyzed the ROI on the dark iris in Fig. 3o of phase 2 because we performed the image processing with Asian patients with dark irises. When the same algorithm was used on a bright-colored iris, the algorithm might have limitation identifying ROI.
We subtracted the RGB average value of iris from the ROI to remove the iris color within the ROI region. However, this approach may not eliminate the noise generated by the iris pattern. To effectively remove iris background, we may obtain the iris data under the ROI separately and subtract the iris value for each pixel (Fig. 2 Phase 4). We used the sclera of the patient's eye for the white color reference. In future analysis, this color variation may be adjusted by using an image or video taken under the same recording conditions using a white paper (i.e., the fiducial marker for color) under the LED light. As an additional method, HSV color space can be adopted to iris color detection phase, which is considered a more effective color space than RGB color space for iris authentication [31]. This approaches would produce more consistent results than using potentially nonidentical sclera color from each patient (Fig. 2 Phase 5).
So far, algorithms have been serialized using only a single core from a quadcore CPU (Intel i7-6700), so it is relatively slow and currently takes approximately 3 min (1.5 s/ frame) to determine the radius and compensate image color with a video clip. We expect to achieve comparable results in less than a minute by using a graphic processing unit (GPU) with parallel processing for executing the algorithm shown in Fig. 2. Color bar values from the theoretical calculation are similar to the patient's LLT values derived by the newly developed algorithm. While only limited number of patients' data were available for our pilot study, the results were reasonable based on the previous study that demonstrated the correlation between the lipid layer thickness and two frequently used dry eye tests by comparing results of tear break-up time and Schirmer's test [32]. Therefore, our approach can be useful in providing objective information to evaluate the lipid layer thickness for patients with dry eye syndrome. For further verification, we are currently analyzing more patients with our algorithm along with their clinical manifestation.

Conclusion
In this study, we utilized white light interference from the lipid layer and implemented a novel algorithm to analyze the thickness of the lipid layer by examining interference color distribution. Our proposed algorithm provides a quantitative measure of the lipid layer thickness even in the presence of eye movement by tracking ROI for each frame. With a small number of patient data set as a feasibility study, our experimental results demonstrate that the lipid layer thickness was different in the subcategories of dry eye patients evaluated by an ophthalmologist. This approach provides a significant step forward in developing a fully automated instrument for evaluation of dry eye category and thus guiding optimal treatment for patients.