Computational approach for tear film assessment based on break-up dynamics

b i o s y s t e m s e n g i n e e r i n g 1 3 8 ( 2 0 1 5 ) 9 0 e1 0 3

Available online at www.sciencedirect.com

ScienceDirect journal homepage: www.elsevier.com/locate/issn/15375110

Special Issue: Innovations in Medicine and Healthcare Research Paper

Computational approach for tear film assessment based on break-up dynamics
ldez b, L. Ramos a,*, N. Barreira a, H. Pena-Verdeal b, M.J. Gira E. Yebra-Pimentel b a b

~a, Spain VARPA Group, Department of Computer Science, University of A Corun Optometry Group, Dept. Applied Physics, Univ. Santiago de Compostela, Spain

article info

Dry eye syndrome is a common disorder of the tear film which affects a remarkable per-

Article history:

centage of the population, impacting on quality of life. The study of the tear film stability is

Received 27 August 2014

essential for the dry eye characterisation. The Break-Up Time (BUT) is a clinical test which

Received in revised form

computes the time the first tear film break-up appears. Besides the time, break-up prop-

8 April 2015

erties can be related to specific aspects of the tear film that could affect dry eye severity.

Accepted 27 April 2015

This work describes a fully automatic methodology to compute the BUT measurement and

Published online 16 May 2015

evaluate the dynamics of break-up areas. This methodology has been tested on a data set consisting of 18 tear film videos, achieving similar results to the manual annotations

Keywords:

marked by the experts. This analysis provides useful additional information for the tear

Tear film

film assessment.

Dry eye syndrome

© 2015 IAgrE. Published by Elsevier Ltd. All rights reserved.

BUT test Video analysis Image processing

1.

Introduction

The tear film is a complex and dynamic trilaminar structure comprising an outer lipid layer, a middle aqueous layer and an inner mucous layer (Holly & Lemp, 1977). It covers the anterior surface of the cornea playing an essential role in the maintenance of ocular integrity health (Guillon, Maissa, & Styles, 1998). Abnormalities in any of the layers can lead to tear dysfunction problems. Concretely, the lipid layer plays a

major role in retarding the evaporation of the tear film during the inter-blink period, and consequently, a deficit of this layer can cause the dry eye syndrome (DES) (Graig & Tomlinson, 2007). The DES is a multifactorial disease of the ocular surface, which affects a significant percentage of the population, and worsens with age Javadi and Feizi (2011); Lemp (2008); Lowther (1977). The prevalence of this syndrome has been increasing in recent years, affecting up to 10e15% of normal population, and 18e30% of contact lenses users. Several factors, such as adverse environmental conditions, use of certain

* Corresponding author. E-mail addresses: [email protected] (L. Ramos), [email protected] (N. Barreira), [email protected] (H. Pena-Verdeal), mjesus.gir
ldez), [email protected] (E. Yebra-Pimentel). [email protected] (M.J. Gira http://dx.doi.org/10.1016/j.biosystemseng.2015.04.009 1537-5110/© 2015 IAgrE. Published by Elsevier Ltd. All rights reserved.

b i o s y s t e m s e n g i n e e r i n g 1 3 8 ( 2 0 1 5 ) 9 0 e1 0 3

medications, or visual tasks that reduce blink rate have contributed to that increment (Fenga et al., 2008; Garcı´aResu´a, Lira, & Yebra-Pimentel, 2005). This disease alters some visual functions such as contrast perception and eye acuity, resulting in symptoms of discomfort and visual disturbance which impact on quality of life (Tutt, Bradley, Begley, & Thibos, 2000; Nichols, Begley, Caffery, & Jones, 1999). The diagnosis of this condition is difficult due to its multifactorial etiology (Khanal, 2008). Normal tear film dynamics require adequate production of tears, retention on the ocular surface, and balanced elimination. Disruption of any of these processes can lead to the condition of dry eye (Tomlinson & Khanal, 2005). Thus, the composition and behaviour of the tear film provide crucial indicators, so the tear film assessment is essential for DES characterisation (Guillon, 1998). There are several clinical tests to evaluate the quality and stability of the tear film on the ocular surface. Among the different tests available, the Break-up Time measure has been widely used in clinical practice. It consists of measuring the time that the tear film remains stable without blinking (Cho et al., 1992; Cho & Brown, 1993; Lee & Kee, 1988). To perform this test, sodium fluorescein is instilled into the eye using a micro-pipette, and the tear film is observed with the help of cobalt-blue filter attached to a slit lamp biomicroscope, and a yellow filter to improve the visibility of the fluorescein emission (Elliott, Fandrich, & Simpson, 1998; Johnson & Murphy, 2005). The patient is instructed to blink three times naturally, without squeezing, in order to distribute the fluorescein over the cornea, and then, he/she maintains the eye open as long as possible (Begley et al., 2006). The BUT is measured as the time elapsed between the last blink and the first appearance of a dark spot on the surface of the cornea, which represents the evaporation of water and the break-up of the tear film. A low BUT measurement corresponds to a limited ocular surface wetting, and it is one of the main signs of an abnormal tear film. The BUT test is affected by low repeatability mainly due to a subjective appreciation of the dark spots, the differences among the experts, and the variability of the tear film. The automation of the break-up assessment would reduce its subjective character, allowing a more accurate evaluation of the tear film stability. A preliminary approach for the automation of the BUT test was conducted in (Cebreiro, Ramos, Mosquera, Barreira, & Penedo, 2011; Ramos et al., 2012). This methodology provides a BUT measurement through an algorithm which consists of locating the different measurement areas from the video sequence, extracting the region of interest, and performing the BUT test over each measurement area. The BUT test only examines the appearance of the first dark spot in the tear film, regardless of the subsequent breakup dynamics. However, the first break-up could appear as a small point or as a large area and its size could increase faster or slower with time. This information is omitted in the BUT test, but it is relevant for understanding tear film instability and its relation to other ocular surface symptoms in dry eye syndrome (Begley et al., 2006). Even though this is an interesting research field, there is no automatic techniques for analysing the break-up evolution in the literature. Therefore, in this work, an automatic methodology is proposed for characterising tear film dynamics over the exposed corneal

91

surface from the emergence of the first break-up in the tear film until the later blink. To this end, the BUT measurement is computed and the break-up areas are segmented in each video frame in order to analyse other break-up parameters, such as size or growing rate. This paper is organised as follows. Section 2 describes the methodology developed for the break-up dynamics assessment. Section 3 summarises the results obtained with the methodology. Finally, Section 4 presents the conclusions and future work.

2.

Methodology

Each tear film video has a duration of several minutes and contains different sequences of interest (SOIs). Each SOI consists of a set of frames between consecutive blinks where the patient maintains the eye open. This way, each SOI contains zero or one BUT measurements. Thus, after the initial full blink which marks the beginning of the SOI, the fluorescein spreads uniformly over the stable tear film. As time passes, the tear film loses stability, and dark areas appear related to the tear film break-up. These rupture zones evolve until the final blink, which marks the end of the SOI. A image processing methodology has been proposed in previous works to automatise this task (Cebreiro et al., 2011; Ramos et al., 2012). This methodology has several stages, as shown in Fig. 1. The first stage consists of extracting the data of interest for break-up evolution assessment (Ramos et al., 2013). It comprises a step for locating the different SOIs and another step for extracting the region of interest (ROI) within each tear film frame. The ROI in each SOI frame corresponds to the visible part of the iris which may vary slightly throughout the sequence depending on the eye aperture and shadows due to outer parts of the eye like eyelids or eyelashes. Once the SOIs are delimited and the ROI is extracted, a preprocessing stage is performed. First, the illumination issues are corrected by means of a contrast normalisation algorithm and, after that, a gray level analysis is performed in order to isolate the break-up areas. Finally, the main stage of this methodology is the tear film dynamics assessment. It consists of measuring the size and growing rate of the breakup areas from the first frame with dark spots until the final blink. For this purpose, the BUT measurement is computed first, and, after that, the break-up areas are segmented and analysed until the end of the SOI. These stages are explained in more detail in the following sections.

2.1.

Selection of sequences of interest

The SOIs forming the tear film video are located as the areas between two consecutive blinks, which delimit the beginning and end of each SOI. These blinks correspond to transitions from closed to open eye and vice versa. If the eye is open, the bright part corresponding to the sclera occupies a significant percentage of the frame, as shown in Fig. 2 (left). Meanwhile, if the eye is closed, the eyelid takes up the entire frame, so, as a consequence, it presents a darker tonality, as shown in Fig. 2 (right).

92

b i o s y s t e m s e n g i n e e r i n g 1 3 8 ( 2 0 1 5 ) 9 0 e1 0 3

Fig. 1 e Stages of break-up dynamics assessment in a tear film video.

Thus, a way to characterise the eye aperture is the use of a measurement which represents the frame tonality. The mean value of intensities Ik represents the average tonality of a frame and can be used to distinguish between open and closed eyes. It is computed as follows: P Ik ¼

i;j Ik ði; jÞ

r*c

(1)

where Ik(i,j) corresponds to the value of the pixel located in the row i and the column j of the frame Ik of the video sequence, and r and c are the number of rows and columns, respectively. Therefore, the detection of blinks is based on detecting intensity variations throughout the sequence. To this end, the finite differences of mean values of gray between consecutive frames are calculated using the following equation: DIk;d ¼ Ikþd  Ik

(2)

where Ik is the mean value of gray for the frame Ik, and d is the distance between the frames for which the difference is computed. This way, DIk;1 represents the differences between consecutive frames whereas DIk;2 stands for the differences between every two frames. Since a blink can last several frames, the computation with a d larger than 1 is required. DIk;d has values close to zero while the eye remains in the same state and higher values when there are variations. Most of the times the blink occurs gradually and the differences between successive frames may not be enough to detect it. In order to detect all the blinks, differences between nonconsecutive frames, separated by a distance of f frames, are f also considered. Thus, DIk is computed as the sum of finite differences between the frame i and the frame i þ d, where d gets values from 1 to f frames:

Fig. 2 e Frames during a blink in a tear film video. Left: Open eye. Right: closed eye.

b i o s y s t e m s e n g i n e e r i n g 1 3 8 ( 2 0 1 5 ) 9 0 e1 0 3

93

Fig. 3 e Automatic detection of SOIs as the intervals between peaks with opposite signs which exceed a length tL. The horizontal red lines represent the adaptive threshold for blink detection. The blue asterisks are the detected blinks used to extract the SOIs marked by dotted blue lines, whereas the black triangles correspond to the start and end of lamp off periods, marked by dotted black lines. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

f

DIk ¼

f X

DIk;l

(3)

l¼1

This sum emphasises the intensity changes produced during a blink. On these differences, a negative peak represents the beginning of a blink, since there is a transition from a lighter frame (open eye) to a darker frame (closed eye). Similarly, at the end of the blink there is a transition from a darker frame to a lighter frame, producing a positive peak, as shown in Fig. 3. In order to identify these peaks, an adaptive threshold tD is defined according to the differences obtained for each video, using the following equation: P 0f nf o DI f tD ¼   k ; I0 k ¼ Ik > 1; k ¼ 1::ðn  dÞ f I0   k

(4)

where DIk;d are the differences of mean values of gray between consecutive frames, n is the number of frames in the tear film video, and d is the distance between the frames where the difference is computed. The range (tD,tD) obtained from this threshold is used to identify transitions from open to closed eye and vice versa. Thus, values outside this range are related to blinks, while values inside this range correspond to the areas where the eye is open. Therefore, the SOIs are identified as those intervals starting with a positive difference and ending with a negative difference. Also, the SOIs should exceed a length tL to ensure that

they have a minimum of frames for break-up assessment. Sometimes there could be two consecutive differences with the same sign. This occurs when the lamp is off or when there are semi-blinks in the SOI. In this case, the lowest absolute values are removed until all pairs of consecutive blinks have opposite signs. Figure 3 shows the areas automatically detected by this approach in a tear film video sequence.

2.2.

Extraction of the region of interest

After locating the SOIs, the next step is extracting the ROI within each frame to discard regions without relevant information. First, the position and size of the iris are identified, since it corresponds to the region where the break-up is evaluated. The eye size does not change through the SOI but its degree of aperture could present slight variations, so different samples of the eye are considered for a robust ROI extraction. Thus, the frames placed at 25%, 50% and 75% are selected for each SOI. The Canny edge detector (Canny, 1986) is applied to the green component of the RGB image of these frames (see Fig. 4, a), and then, the dark-to-light edges are discarded, since they are related to the sclera area (see Fig. 4, b). Assuming the iris has an approximately circular shape, the segmentation is carried out by correlation in the frequency domain between the edge frames and a set of circular and elliptical masks with different radii covering typical eye sizes. Therefore, the radii and position of the mask with the highest matching value delimit the iris, as Fig. 4, c.

94

b i o s y s t e m s e n g i n e e r i n g 1 3 8 ( 2 0 1 5 ) 9 0 e1 0 3

Fig. 4 e ROI extracted by correlation. (a) Edge image obtained by Canny edge detector (b) Edge image discarding incoming pixels (c) Eye size obtained from the best correlation value.

Moreover, the eye presents slight motions throughout the video, so it is necessary to register the ROI in the video sequence. An alignment by correlation is performed in the frequency domain between the edge images of each frame and the immediately preceding frame. After this, the mask selected in the previous step is applied to the aligned frames. This way, the methodology is independent of slight motions of the ROI, since the edge alignment between consecutive frames produces a good enough match. In some cases, the eye is not fully open and the ROI contains outer parts like eyelids or eyelashes. These regions do not contain relevant information for the analysis and could mislead the results. For these reasons, a further adjustment is performed in each SOI frame, as shown Fig. 5. It consists of cropping the top and the bottom of the circular ROI, according to the features of each eye. For this purpose, the number of edges is computed in each row i, as calculation in Eq. (5). Then, the upper limit is selected from the furthest row in the upper half of the image that accumulates more edge points than a variable threshold, as calculation in Eq. (6). In the same way, the lower limit is the closest row in the lower half of the image

that accumulates more edge points than a variable threshold, as calculation in Eq. (7). AccðiÞ ¼

edgesði; jÞ

(5)

j¼i

  lupper ¼ max i; AccðiÞ > tupper

(6)

llower ¼ minfi; AccðiÞ > tlower g

(7)

The thresholds are computed as a percentage of the maximum number of edge points found in each half of the image. Given the nature of the eye, the upper eyelid and eyelashes usually invade more ROI than the lower, so the parameters a and b are used as weights for the upper and lower threshold, as calculation in Eqs (8) and (9). tupper ¼ amaxfAccðiÞ; 0  i < rows=2g

(8)

tlower ¼ bmaxfAccðiÞ; rows=2  i < rowsg

(9)

Furthermore, the radius is slightly reduced to get rid of noise at the boundaries of the iris.

2.3.

Fig. 5 e Adjustment of the ROI to discard irrelevant information. The top and the bottom are cropped, and the radii are slightly reduced.

cols X

Frame preprocessing

The tear film images are non-uniformly illuminated, so there are luminosity and contrast heterogeneity within the images. This problem could affect the break-up characterisation, mistaking poorly illuminated areas for real break-up areas. In order to overcome this problem, a lighting correction based on the approach proposed by Foracchia et al. (Foracchia, Grisan, & Ruggeri, 2005) is performed. This process consists of normalising the lightness and contrast variability in images, based on estimating both features in background small areas, spreading to the whole image and then removing from it. Figure 6 shows an example of this step, where the preprocessed frame presents uniform luminosity. After the luminosity correction, the break-up is detected from the intensities of the green component of the normalised tear film frames. Each tear film video presents different illumination and amount of fluorescein instilled, so not all the SOIs have the same intensity levels. Furthermore, the dark pixels at the break-up vary in a range close to zero according to lighting conditions, but not exactly zero. For these reasons, a threshold is set to determine the break-up threshold tb in

b i o s y s t e m s e n g i n e e r i n g 1 3 8 ( 2 0 1 5 ) 9 0 e1 0 3

95

Fig. 6 e Normalisation by contrast and luminosity correction. Left: original frame. Right: preprocessed frame where the lighting is uniform due to the normalisation.

each SOI, that is, the maximum intensity of a pixel to be considered as part of the break-up. For this purpose, two alternatives have been explored. In the first approach, tb is computed from the cumulative histogram of a reference frame located at the beginning of the SOI. The cumulative histogram of this reference frame is built and the gray levels of a percentage pb of the darkest pixels are analysed, as shown in Fig. 7. The break-up threshold tb corresponds to the largest value of these gray levels. Pixels with values below the break-up threshold are considered as part of the break-up area. Thus, the break-up areas can be segmented, as shown in Fig. 8. The second approach uses the frame which delimits the end of the SOI. This frame allows a better characterisation of the break-up pixels because it represents the maximum

expansion of the break-up in the SOI. In this frame, a multilevel thresholding (Arora, Acharya, Verma, & Panigrahi, 2008) from the statistical distribution of the intensity values is applied. This method uses the mean and the variance of the image to find optimum thresholds for segmenting the image into multiple levels. The algorithm 1 summarises the steps for multilevel thresholding. It takes initially the range which include all the intensities, that is, from 0 to 255 and computes the mean m and the standard deviation s of all the pixels in this range. From these values, the limits t1 and t2 for a new sub-range are calculated as follows: t1 ¼ m  ks t2 ¼ m þ ks

(10)

where k is a free parameter. The algorithm is applied recursively on sub-ranges computed from the previous step so as to

Fig. 7 e Cumulative histogram of the reference frame. The break-up threshold tb is obtained as the highest gray level of a percentage pb of the darkest pixels.

96

b i o s y s t e m s e n g i n e e r i n g 1 3 8 ( 2 0 1 5 ) 9 0 e1 0 3

Fig. 8 e Break-up segmentation. Left: original frame at the end of a SOI. Right: break-up areas detected by thresholding using tb. find a threshold and a new sub-range for the next step, until the number of chosen thresholds is reached. Algorithm 1. Multilevel thresholding.

a better definition of the input frame, but after reaching a certain number of levels, the addition of classes hardly varies the final segmentation. The threshold tb is computed from lb, one of the lowest levels obtained in the multilevel thresholding, since they represent the darkest areas of the frame. Thus, these levels are analysed in order to get the best threshold to segment the break-up areas. The threshold tb corresponds to the upper limit of the selected level lb. Figure 10 shows the break-up segmentation by selecting different levels lb as threshold tb. The lower levels segment the darkest areas whereas the upper levels segment almost all the image.

2.4.

Figure 9 shows the frame segmentation using different thresholding levels. Increasing the number of classes provides

Break-up dynamics

The main stage of the methodology is the break-up dynamics assessment, which consists of analysing break-up evolution from the BUT measurement until the following blink. Thus, in the first step, the BUT is computed as the time elapsed since the last blink until the tear film break-up, that is, the first frame where dark areas appear on the surface of the eye. In order to detect the emergence of these points, the pixels with intensities lower than tb are counted in each frame. When the amount of these pixels exceeds a threshold te, this frame is considered as a break-up frame.

Fig. 9 e Multilevel thresholding with different levels. From left to right, top to bottom: Original frame, 4-level, 6-level, 8-level, 10-level, and 12-level thresholding.

b i o s y s t e m s e n g i n e e r i n g 1 3 8 ( 2 0 1 5 ) 9 0 e1 0 3

97

Fig. 10 e Break-up segmentation by multilevel thresholding. The top left frame is the original result of multilevel thresholding and the rest are the break-up segmentation with tb obtained with values of lb from 1 to 5.

An evolution curve is built with the percentage of these pixels from the eye fully open until the following blink. Small variations produce curves with irregular slopes, so a second order polynomial is fitted to discard these fluctuations, as shown in Fig. 11. In some cases, the gradient of the curve is zero because the tear does not break in the interval, as shown in Fig. 11 (top). By contrast, if there is a break-up, the percentage of black increases with the time since the fluorescein is not regenerated, as shown in Fig. 11 (bottom). The threshold te is obtained from a percentage pe of the total height of the evolution curve. The BUT measurement is computed as the time elapsed from the beginning of the SOI until the curve exceeds this threshold, as shown in Fig. 11 (bottom) (Ramos et al., 2014). Once the BUT measurement is computed, the break-up evolution is examined from the break-up frame until the end of the SOI. For this purpose, nd frames located every 0.5 s (15 frames for a 30 fps video) after the first break-up until the last blink are selected. For each of these frames, the area BUk of the tear film break-up is obtained from the percentage of break-up pixels in the ROI, as follows: BUk ¼

nkbu 100 nROI

(11)

where nkbu is the break-up pixels of the frame Ik and nROI is the total number of pixels of the ROI, which is the same for all the frames in the SOI. The growing rate rb is computed from the first break-up until the final blink by calculating the mean of the differences between the break-up areas segmented in the selected frames, as follows: rb ¼

Pnd

k¼1 DBUk

nd

DBUk ¼ BUkþ15  BUk

(12)

(13)

This analysis provides additional information to the BUT test since the break-up size and its evolution are quantified from the first appearance until the later blink. Therefore, this

assessment allows characterisation of whether the break-up is small or large and whether it evolves quickly or slowly.

3.

Results

The methodology has been tested on a dataset of tear film
videos provided by the Facultade de Optica e Optometrı´a of the Universidade de Santiago de Compostela. These videos have been recorded with a Topcon DV-3 camera attached to the Topcon SL-D4 slit lamp. The dataset consists of 18 videos from healthy patients with ages ranging from 19 to 33, varying from very dry eye to an eye with no visible dryness. Patients were requested to keep the eye open until they felt a real need to blink. This video dataset contains 113 SOIs which have been manually annotated by four experts. Furthermore, a selection of 140 random frames in different points of break-up dynamics, has been manually segmented by two different experts. First, the SOI location is validated by comparing the SOIs detected by the system with the SOIs delimited by the experts. Table 1 shows the performance of SOI location obtained with values of f from 1 to 3. As can be seen, using f ¼ 1 some SOIs are missed as well as some invalid sequences are detected, whereas including differences between non consecutive blinks provides better results, since this sum emphasises the intensity changes produced during a blink. This way, with f ¼ 2 nearly all the SOIs are located and almost no invalid sequence is included. However, with f  2 the improvement is negligible, so f is set to 2 since the computational complexity is lower and this value works well detecting SOIs with BUT measurement as well as discarding invalid sequences which contains semi-blinks or periods where the lamp is off. In order to validate the break-up assessment, the parameters are adjusted by building ROC curves and selecting the best values on the basis of sensitivity and specificity criteria (Fawcett, 2006; Powers, 2007). Due to the limited size of the

98

b i o s y s t e m s e n g i n e e r i n g 1 3 8 ( 2 0 1 5 ) 9 0 e1 0 3

Fig. 11 e Evolution curves (dotted green line) fitted by a second order polynomial function (continuous red line). Top: flat curve due to a stable tear film. Bottom: the black area increases with time and the BUT is detected when the curve exceeds the threshold te. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1 e Performance of SOI detection with different values of the parameter f. f¼1 f¼2 f¼3

% located SOIs

Invalid sequences

92.32 97.14 97.14

7 (semi-blinks, lamp off) 1 (semi-blink) 1 (semi-blink)

dataset, 10-fold cross-validation is used in the experiments to assess the generalisation capability (Rodriguez, Perez, & Lozano, 2010). Therefore, the dataset is divided into 10 parts and an iterative process is carried out 10 times. At each iteration, one of the parts is used as a test set, and the remaining parts are used as a training set. Finally, ROC curves are built for each 10 k-fold training set, and analysed to adjust the value for each parameter. Once each parameter is set to the best value, the accuracy is extracted from the k-fold test sets.

b i o s y s t e m s e n g i n e e r i n g 1 3 8 ( 2 0 1 5 ) 9 0 e1 0 3

Fig. 12 e Break-up areas manually annotated by two experts in the same frame. Red regions are related to the areas marked by Expert 1 whereas blue regions correspond to Expert 2 annotations. The intersection between the regions segmented by the both experts is represented in green.

First, the accuracy of the break-up segmentation approaches is tested. To this end, each pixel is classified as break-up or background, so that values for sensitivity and specificity are obtained by comparing the matches pixel by pixel between the manual annotations and the automatic results. Note that labelling the break-up areas by hand is a tedious and subjective task so that there is a remarkable

99

disagreement between the experts, as shown in Fig. 12. The coincidence of the delimited break-up areas is around 60%, so the intersection between both experts is used for adjusting the parameters. Two alternatives have been proposed for computing the break-up threshold. In the first approach, the break-up threshold tb is obtained from the percentage pb of the darkest pixels in the cumulative histogram of the reference frame. A ROC curve is built with pb taking values from 5% to 95% for each k-fold training set. The best parameter would yield a point in the upper left corner of the ROC space, representing 100% sensitivity and 100% specificity. In order to get a compromise between these two measures, a line from the upper left point to the opposite corner is drawn, and the pb is set to the value of the ROC curve crossed by this line. The best values for pb are always between 40% and 45% in the different ROC curves corresponding to the 10-fold training sets, so a mean ROC curve is computed as a representation of the training step. This ROC curve is built from the mean values of specificity and sensitivity in the 10 k-fold training sets, and used to set the value of pb, as shown in Fig. 13. Thus, pb is set to 45% of the darkest pixels of the cumulative histogram. This value was checked in the 10 k-fold test sets providing a mean accuracy of 82.71%, from 73.26% in the worst case to 89.46% in the best case. The second approach uses the frame which delimits the end of the SOI. The threshold tb is computed from upper limit of lb, which is one of the levels obtained by the multilevel thresholding. In order to select the value of lb which provides the best segmentation, all the levels of 4level, 6level, 8level, and 10level thresholding have been analysed for the different k-fold training sets. The different 10 k-fold training sets provide the same results for the best lb in each thresholding.

Fig. 13 e Mean ROC curve for the 10 k-fold training sets obtained with percentages from 5% to 95% of the darkest pixels. Sensitivity and specificity were computed by comparing the matches between the automatic segmentation and the intersection of the manual annotations performed by both experts.

100

b i o s y s t e m s e n g i n e e r i n g 1 3 8 ( 2 0 1 5 ) 9 0 e1 0 3

Fig. 14 e Mean ROC curve for the 10 k-fold training sets obtained with values of lb from the levels of 4¡level, 6¡level, 8¡level, and 10¡level thresholding. Sensitivity and specificity were computed by comparing the matches between the automatic segmentation and the intersection of the manual annotations performed by both experts.

Thus, mean ROC curves computed from the mean values of sensitivity and specificity in the 10 k-fold training sets are used as a representation of the training iteration, as shown in Fig. 14. The 6level presents results slightly higher than the

4level, 8level, and 10level thresholding. Therefore, lb is set to level 3 in the 6level thresholding, since it provides the best compromise between sensitivity and specificity. This value was checked in the 10 k-fold test sets providing a mean

Fig. 15 e Mean ROC curve for the 10 k-fold training sets obtained with percentages from 5% to 95% of the total height of the evolution curve. Sensitivity and specificity were computed by comparing the matches between the expert average and the automatic BUT measurements.

b i o s y s t e m s e n g i n e e r i n g 1 3 8 ( 2 0 1 5 ) 9 0 e1 0 3

101

Fig. 16 e Break-up dynamics in two different SOIs. Top: the first break-up represents the 33% of the whole ROI and evolves with a growing rate of 1.75% per second until reach the 47% of the whole ROI. Bottom: the first break-up covers the 11% and the growing rate is about 7.5% per second, occupying almost the 70% of the ROI at the end of the SOI.

accuracy of 83.78%, from 79.55% in the worst case to 89.97% in the best case. These approaches get acceptable results considering the high variability between the experts. The sensitivity and specificity provided by the multilevel thresholding is slightly better, so tb is computed from the upper limit for the level lb 3 in the 6level thresholding. This break-up threshold tb will be used in the step for computing the BUT measurement.

In the same way, the parameter pe related to the BUT measurement is set by comparing the break-up frames automatically classified with respect to the average of the values provided by 4 different experts. Thus, frames located before the BUT are considered as frames without break-up and the subsequent frames after the BUT are considered as frames with break-up. A ROC curve is built with pe taking values from 5% to 95% of the total height of the evolution curve for each k-

102

b i o s y s t e m s e n g i n e e r i n g 1 3 8 ( 2 0 1 5 ) 9 0 e1 0 3

fold training set. The best values for pe are always between 50% and 55% in the different k-fold iterations, so the ROC curve computed from the mean values of specificity and sensitivity in the 10 k-fold training sets is used to set the value of pe, as shown in Fig. 15. Thus, pe is set to 55%, value which provides a mean accuracy of 79.68%, being 76.16% the worst case, and 84.88% the best case in the 10 k-fold test sets. Once the parameters have been adjusted, the break-up dynamics is evaluated. For this purpose, the first break-up frame and the following frames every 0.5 s until the end of the SOI are selected. The percentage of the break-up pixels over the whole ROI is computed for each of these frames. Furthermore, the growing rate is computed from the mean of the differences among the break-up areas in the selected frames, second by second. This analysis provides additional information to the BUT test since the break-up size and its evolution are quantified from the first appearance until the later blink. Therefore, this assessment allows characterisation of whether the break-up is small or large and whether it evolves quickly or slowly. Figure 16 shows two examples where the break-up dynamics graphs distinguish between a large break-up with slow and small increase, and a small break-up with massive expansion.

4.

Conclusion and future work

This work describes a fully automated methodology for characterising tear film dynamics over the exposed corneal surface from the BUT measurement until the later blink. The procedure first extracts the data of interest for break-up evolution assessment by locating the different SOIs and defining the ROI within each SOI frame. Then, a preprocessing stage is performed to normalise the contrast and illumination of the tear film frames and determine a suitable threshold to segment the break-up areas. Finally, the tear film dynamics are analysed. In this stage, the break-up time is computed and the break-up areas are characterised in terms of size and growing rate. This methodology has been tested on a data set consisting of 18 tear film videos with 60 SOIs achieving similar results to the manual annotations marked by the experts. The characterisation of the break-up dynamics allows a quantitative, objective analysis of tear instability, as an extension of BUT measurement, which focuses only on the time. Besides the subjectivity, the characterisation by hand of break-up dynamics is hard to obtain so this analysis provides useful additional information to clinical practice. Future work in this field includes an analysis of the repeatability of the BUT test and break-up dynamics between different clinical revisions.

Acknowledgements This research has been economically supported in part by the
n e Ordenacio
n Universitaria Consellerı´a de Cultura, Educacio of the Xunta de Galicia through the agreement for the Singular Research Center CITIC.

references

Arora, S., Acharya, J., Verma, A., & Panigrahi, P. K. (2008). Multilevel thresholding for image segmentation through a fast statistical recursive algorithm. Pattern Recognition Letters, 29(2), 119e125. Begley, C., Himebaugh, N., Renner, D., Liu, H., Chalmers, R., Simpson, T., et al. (2006). Tear breakup dynamics: a technique for quantifying tear film instability. Optometry & Vision Science, 83, 15e21. Canny, J. (1986). A computational approach to edge detection. IEEE PAMI, 679e698. Cebreiro, E., Ramos, L., Mosquera, A., Barreira, N., & Penedo, M. (2011). Automation of the tear film break-up time test. In ISABEL. Cho, P., & Brown, B. (1993). Review of the tear break-up time and a closer look at the tear break-up time of Hong Kong Chinese. Optometry & Vision Science, 70, 30e38. Cho, P., et al. (1992). Reliability of the tear break-up time technique of assessing tear stability and the locations of tear break-up in Hong Kong Chinese. Optometry & Vision Science, 69, 879e885. Elliott, M., Fandrich, H., & Simpson, T. (1998). Analysis of the repeatability of tear break-up time measurement techniques on asymptomatic subjects before, during and after contact lens wear. Contact Lens & Anterior Eye, 21, 98e103. Fawcett, T. (Jun. 2006). An introduction to ROC analysis. Pattern Recognition Letters, 27(8), 861e874. Fenga, C., Aragona, P., Cacciola, A., Spinela, R., Nola, C. D., Ferreri, F., et al. (2008). Melbonian gland dysfunction and ocular discomform in video display terminal workers. Eye, 22, 91e95. Foracchia, M., Grisan, E., & Ruggeri, A. (2005). Luminosity and contrast normalization in retinal images. Medical Image Analysis, 9(3), 179e190.
n Garcı´a-Resu´a, C., Lira, M., & Yebra-Pimentel, E. (2005). Evaluacio
venes universitarios. Revista Espan ~ ola de superficial en jo Contactologı´a, 12, 37e41. Graig, J., & Tomlinson, A. (2007). Importance of the lipid layer in human tear film stability and evaporation. Optometry & Vision Science, 74, 8e13. Guillon, J. (1998). Non-invasive tearscope plus routine for contact lens fitting. Contact Lens & Anterior Eye, 2, S31eS40. Guillon, M., Maissa, C., & Styles, E. (1998). Relationship between pre-ocular tear film structure and stability. Lacrimal Gland, Tear Film, and Dry Eye Syndromes, 438, 401e405. Holly, F., & Lemp, M. (1977). Tear film physiology and dry eyes. Survey of Ophthalmology, 22, 69e87. Javadi, M.-A., & Feizi, S. (2011). Dry eye syndrome. Journal of Ophthalmic & Vision Research, 6(3). Johnson, M., & Murphy, P. (2005). The effect of instilled fluorescein solution volume on the values and repeatability of TBUT measurements. Cornea, 24, 811e817. Khanal, S. (2008). Dry eye diagnosis. Investigative Opthalmology & Visual Science, 49, 1907e1914. Lee, M. J. H., & Kee, M. C. W. (1988). The significance of tear film break-up time in the diagnosis of dry eye syndrome. Korean Journal of Opthalmology, 21, 69e71. Lemp, M. (2008). Advances in understanding and managing dry eye disease. American Journal of Ophthalmology, 46, 350e356. Lowther, G. (1977). Dryness, tears and contact lens wear. USA: Reed Elsevier. Nichols, K., Begley, C., Caffery, B., & Jones, L. (1999). Symptoms of ocular irritation in patients diagnosed with dry eye. Optometry & Vision Science, 76, 838e844. Powers, D. M. W. (2007). Evaluation: from precision, recall and Ffactor to ROC, informedness, markedness & correlation. Journal of Machine Learning Technologies, 2, 37e63.

b i o s y s t e m s e n g i n e e r i n g 1 3 8 ( 2 0 1 5 ) 9 0 e1 0 3


s, M., PenaRamos, L., Barreira, N., Mosquera, A., Curra
dez, M., et al. (2013). Computational approach Verdeal, H., Gira for measuring the tear film break-up time in an unsupervised manner. Advanced Techniques for Knowledge Engineering and Innovative Applications, 246, 254e267.
s, M., PenaRamos, L., Barreira, N., Mosquera, A., Curra
ldez, M., et al. (2012). Adaptive parameter Verdeal, H., Gira computation for the automatic measure of the tear break-up time. In 16th International Conference on Knowledge-Based and Intelligent Information & Engineering Systems (KES'12), 243 pp. 1370e1379). Ramos, L., Barreira, N., Mosquera, A., Penedo, M. G., YebraPimentel, E., & Garcı´a-Resu´a, C. (2014). Analysis of parameters

103

for the automatic computation of the tear film break-up time test based on cclru standards. Computer Methods and Programs Biomedicine, 113(3), 715e724. Rodriguez, J., Perez, A., & Lozano, J. (2010). Sensitivity analysis of k-fold cross validation in prediction error estimation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(3), 569e575. Tomlinson, A., & Khanal, S. (2005). Assessment of tear film dynamics: quantification approach. Ocular Surface, 3(2), 81e95. Tutt, R., Bradley, A., Begley, C., & Thibos, L. (2000). Optical and visual impact of tear break-up in human eyes. Investigative Ophthalmology & Visual Science, 41, 4117e4123.