Automatic dynamic tear meniscus measurement in optical coherence tomography

An image processing algorithm is developed for quantitative assessment of tear meniscus dynamics from continuous optical coherence tomography (OCT) measurements. Clinical utility of dynamic OCT tear meniscus measurement is assessed in studies of tear meniscus parameters. The results indicate that any apparent changes in the early post-blink phase meniscus parameters are essentially related to the longitudinal movements of the eye and not to the formation of tear meniscus corresponding to tear film build-up. Dynamic acquisition of tear film meniscus is essential for providing reliable estimates of its parameters such as height, depth, and area.


Introduction
Tear film is important to the health and optics of the human eye. Measurement of tear film properties provides clinicians with investigative capabilities in the evaluation of ocular surface health [1,2], dry eye diagnosis [3][4][5], or the assessment of contact lens wear [6,7].
Depending on the approach taken, tear film measurement techniques can be divided into invasive and noninvasive [8,9], subjective and objective [10], and into static and dynamic [11]. Measurement of tear meniscus parameters such as its height, depth, area, and calculation of its volume can be used for assisting dry eye diagnosis [12], assessing the efficacy of artificial tear supplements [13], evaluating the effects of contact lens wear [14,15], and facilitating, using mathematical models, estimation of tear film thickness [16,17]. Methods of tear meniscus measurement include en-face slit lamp observation, reflective meniscometry [18][19][20] and cross-sectional measurement with optical coherence tomography (OCT) [21][22][23][24].
The latter belongs to the class of noninvasive methods concerning tear meniscus parameters that has the potential to be conducted objectively and in a dynamic fashion.
The majority of methods for extracting tear meniscus parameters from OCT images require manual input from an operator [23,25]. To the best of our knowledge the only method that was able to perform the task of extracting tear meniscus parameters in an automatic fashion was proposed by Tsechpenakis and Wang [26]. The method, based on conditional random field theory, appears to be accurate but the result presented in that work was for a static single measurement and there was no indication how the method performs for a set of dynamic OCT images acquired immediately after a blink. As the authors conclude, an automatic quantitative estimation of the tear menisci is a challenging task.
The repeatability of OCT based tear meniscus measurement has been considered but it was hindered by the necessity of the manual input from an operator. Johnson and Murphy [27] studied the tear meniscus height using two static OCT measurements on two different days. Wang and colleagues [28] used dynamic acquisition of tear meniscus at eight frames per second. However, to evaluate the repeatability of tear meniscus parameters, only static measurements (the first good image in each acquisition) were taken into consideration. Bitton and colleagues [29] assessed the variability of tear meniscus height using 10 consecutive scans using two graders but they had to reevaluate their results after a training session among the observers. Hence, manual estimation of the tear meniscus parameters is subjective and very laborious.
The dynamics of OCT based tear meniscus measurement has not been fully investigated particularly in the sense of tear meniscus formation as well as the reliability of a single postblink static acquisition. Measurement of tear meniscus dynamics can be conducted with commercially available instrumentation. Many OCT systems have the capability of performing multiple B-scan measurements at the standard video rate of 24 frames per second. At that sampling frequency, an average blink of 250 ms is confounded to about six frames [30] while the tear film build-up phase interlocked with further slower post-blink movement of the upper eyelid may take additional 12 frames [31,32].
Until now, there were no automatic techniques that could facilitate quantitative assessment of tear meniscus dynamics with OCT. In this paper, we develop an automatic image processing procedure for dynamic OCT based tear meniscus measurement and assess its clinical utility in studies of tear meniscus parameters.

Measurement protocol
B-scans of the inferior tear meniscus (see Fig. 1) were obtained using commercially available Spectral Optical Coherence Tomography instrument -SOCT Copernicus HR (Optopol, Poland). The instrument offers several scanning protocols for the anterior eye segment, from which the animation mode was chosen. This protocol allows observing the dynamics of tear meniscus. Setting examination parameters include the scanning width, number of A-scans per B-scan, the total number of B-scans and the scanning angle. The scanning angle was set to 90 degrees. As the region of interest is relatively small, the scanning width was set to the minimum available value of 4 mm. This parameter does not affect other parameters. The maximum number of A-scans per B-scan in this mode is 1800 and it determines the highest possible image quality. The maximum possible number of B-scans, up to which the number of A-scans did not decrease, was set to 90.
Ten volunteers including six females and four males, aged from 24 to 51 y.o., have been recruited for the study from the staff and students of Wroclaw University of Technology. All subjects gave informed consent and the study adhered to the tenets of the Declaration of Helsinki. There were no exclusion criteria as the purpose of the measurements was to obtain a range of different physiological characteristic of tear meniscus dynamics. Each subject was asked to place his chin on the chin rest and lean his forehead on the headrest. To measure the inferior tear meniscus, the subject was asked to look straight ahead (with the help of the fellow eye and additionally placed mirrors on the sides of the instrument) and refrain from any head or eye movements. The orientation of the B-scan plane was central (with respect to the iris outline) to the posterior region of the eye-eyelid junction and normal to the eyelid.
Also, to capture the tear meniscus formation subject was asked to blink once, right after the scanners began to work. If the subject moved the eye or blinked too early or too late, the procedure was repeated. Every subject was measured at least 10 times. No additional eye drops were used. All measurements were conducted in the same environmental conditions. Room temperature and humidity were measured with thermo-hygrometer (C2131, Comet, Czech Republic). The room temperature was (mean ± standard deviation) 23.5 ± 0.3°C while the humidity was 35 ± 3%. The measurements were performed in mesopic conditions. The sampling frequency of the acquired B-scans was approximately 22.7 Hz. In every single measurement 90 B-scans were obtained. They were then processed off-line by the proposed custom written algorithm in a backward loop beginning from the last scan.

Image processing procedure
The SOCT Copernicus HR is primarily designed for posterior segment imaging. The ability to image an anterior segment is achieved with an additional corrective lens module that is attached to the head of the instrument. The acquired images of the anterior segment are therefore distorted. Our aim was to design an image processing methodology for extracting information on tear meniscus from the set of dynamically acquired anterior segment images that is general and not restricted to a particular commercially available instrument. Hence, all the parameters of the tear meniscus estimated through the procedure will be given in pixels rather than in their corresponding physical dimensions.
The size of single B-scan is 900x1009 pixels. The region of interest (ROI) is extracted from the original image, in every iteration starting form the last scan, and it has an empirically chosen size of 351 × 201 pixels. The normalized correlation function is used to track the ROI position from a scan to a scan. Since the OCT images contain substantial amount of noise in both background and imaged structures that change from an image to image, applying of-theshelf image processing techniques did not result in a successful estimation of the tear meniscus parameters.
The proposed tear meniscus parameter estimation algorithm is depicted in Fig. 2. Using an image with ROI, the first step of the algorithm is generation of a binary mask that separates unwanted interference and meniscus profile information from the image components containing eye and eyelid structure. This mask is then used in the algorithm that estimates the Eye-Eyelid Edge (EEE) while its complement is used for estimating the Tear Meniscus Profile (TMP). Information on EEE and TMP is then merged to produce estimates of tear meniscus area (TMA), tear meniscus depth (TMD), and tear meniscus height (TMH). Mask generation is realized through three important steps (see Fig. 3): filtering, multithreshold segmentation, and class merging. Filtering is performed by convolving the original image with a circular mask of empirically chosen radius of two. Multi-threshold segmentation is performed with Otsu method for five classes [33]. Class merging transforms the five level output of the Otsu's algorithm, pre-filtered with a median filter that ensures discontinuities in the meniscus profile, to a binary image, in which the last four classes are combined. Finally, an addition of a U-frame and the flood-fill operation is performed. This particular procedure of the binary mask was dictated by the statistics and the structural content of the ROI image, where the same intensities and similar structures were present in both unwanted and useful signal areas. The smallest number of classes in the multi-thresholding segmentation that would successfully separate the two areas was sought and empirically set to five. Most of the unwanted interference was then confound to class 1, while the useful signal was distributed among classes 2 to 5. Hence, merging of the latter resulted in a binary mask that was further employed for EEE and TMP extraction. The binary mask was directly used to extract EEE by column-wise searching for the last zero-value pixel corresponding to the estimate of the edge at that location. The edge pixels were divided into two parts representing the eye edge and the eyelid edge. To provide a smooth estimate of the EEE, both parts were separately fitted by the fifth and the second order polynomial expansions, respectively. Finally they were both merged again to create the EEE (see Fig. 4). To extract TMP information, the ROI image is first multiplied by an eroded complement of the binary mask. The erosion was performed to rule out high intensity pixels (i.e., those with the intensity comparable to that of TMP pixels) possibly protruding from the eye and eyelid structure to the mask area. Further, similarly to the steps employed for the mask generation itself, average filtering is performed on the image in order to strengthen useful information and suppress the unwanted interference. This is followed by Otsu segmentation to separate TMP pixels from the noise. Since the intensity statistics in this area of ROI substantially change from measurement to measurement, the number of classes in Otsu method needs to be selected conditionally. Given the particular OCT images used in study, the area in which the possible TMP pixels could be found is estimated to contain no more than 200 pixels. Thus if the number of pixels in the highest Otsu class is greater than 200, the number of classes is increased by one and the condition is checked again. Empirically it was observed that the class number properly extracting the TMP information was never less than three. Figure 5 shows a representative example of extracted TMP information (denoted by gray dots). The TMP information, in terms of pixels is then inspected for outliers. Further, for each horizontal pixel value, a median of the corresponding vertical pixel values was calculated. This ensured a unique x to y characterization of the TMP that could be parametrically modeled, for example with an even order polynomial expansion. Since many measurements do not contain meniscus information at the junctions with eye and eyelid surfaces, and to avoid extrapolation, additional boundary points from the EEE are taken before TMP fitting is performed. Those boundary points were set to fulfill the condition of Y EEE < min(Y TMP ). It is worth noting that a polynomial model is arbitrary. It is likely that the water surface of the TMP at the contact with the tissue would be tangent. Hence, a circle tangent to the EEE would probably be the model of choice in case the images were appropriately corrected for distortions.
The final step of the image processing procedure consists of merging the fitted EEE and TMP curves and estimating the position of meniscus-eyelid junction (point A, see Fig. 2), meniscus-eye junction (point C), and the eye-eyelid junction (point B). Additionally, point D is estimated using a bisector of the eye-eyelid angle. Estimated positions of points A, B, C, and D are then used to estimate the following tear meniscus parameters: • TMA -defined as a number of pixels contained within the ABC area enclosed by the EEE and TMP, • TMD -defined as the distance between points B and D, • TMH -defined as the distance between points A and C. The image processing procedure is repeated for all post-blink scans and results in timevarying estimates of TMA, TMD, and TMH.

Results with discussion
In general, time-varying estimates of TMA, TMD, and TMH revealed little variation in the post-blink meniscus dynamics. An example in which the post-blink increase in tear meniscus parameters is clearly visible is shown in Fig. 6. In particular, it was difficult to regularly capture the phase of tear meniscus formation. Nevertheless, the average values of TMA, TMD, and TMH were calculated over about two and half second period in which the first 0.5 second was omitted to account for any tear film build-up and post-blink longitudinal movement of the eye. This corresponded to 50 scans that were taken from each series into consideration.
The subject average (over multiple measurements) of TMA, TMD, and TMH are shown in Fig. 7. To keep all the results within the same scale and units of pixels, square root of TMA is shown. Error bars denote one standard deviation. The result of the Levene's variance test (p < 0.01) prohibited the use of standard ANOVA. The Kruskal Wallis revealed statistically significant differences (p<0.001) between the subjects, indicating that the OCT based tear meniscus measurement is able to relatively precisely characterize tear meniscus parameters of an individual.  To ascertain the source of variability (or lack of) in the post-blink phase of tear film meniscus, an experienced subject was asked to exercise regular tiny eye twitching during the measurement. Twitching was evident in slight movement of the upper and lower eyelid and corresponded to subject-induced longitudinal eye movements of amplitude much great then than naturally occurring in the eye during steady viewing conditions [34]. An example of the results from this exercise is depicted in Fig. 8. Note that the changes in magnitude of TMA, TMD, and TMH are substantial. However, similar variability but less regular were observed in the earlier standard measurements of tear meniscus parameters. The proposed algorithm for automatic extraction of tear meniscus parameters from OCT images has some limitations. Figure 9 shows a collection of different types of OCT tear meniscus measurement acquisitions encountered in the study. In the top row there are two images with complete TMP information, one with discontinued TMP, and two with partial TMP information. Ticks and crosses in each image indicate whether the proposed algorithm could or could not successfully process the particular type of acquisition. Placing both signs in an image indicate a probabilistic behavior of the algorithm resulting in a noisy time series of TMA, TMD, and TMH. In the bottom row of Fig. 9 there are three images with specular reflex (weak and central, weak with small meniscus, and strong and non-central), an image with tear debris, and an image of a subject with irregular scleral topography. The proposed procedure of automatic meniscus extraction could successfully resolve the majority of images with specular reflexes unless it was very strong.

Conclusions
The parameters of tear meniscus depend, aside from the physical parameters of the tear film itself, on several external factors such as the temperature and humidity [35] of the environment as well as subject related factors such as the form of the blink [36], palpebral aperture, and the morphology of the cornea and eyelid junction [37]. The environmental factors can be relatively easily maintained. However, the subject related aspects are difficult to regulate.
OCT based B-scan tear meniscus measurement is important for the assessment of tear parameters because the volume of tear film can be related to the resulting tear meniscus parameters [38]. It can help in the diagnosis of dry eye and the assessment of contact lens wear.
A methodology for automatic quantitative extraction of tear meniscus parameters has been proposed. The method performs well for different types of eye-eyelid junction morphology, a wide range of post blink dynamics, and a varied size of tear meniscus. The method is also robust to moderate amounts of interference resulting from the specular reflex.
The results of this study can be easily generalized to the measurements of upper tear meniscus parameters. The procedure of image processing is equally applicable with the difference of flipping the acquired image of the upper tear meniscus about the central vertical axis.
The results of the study indicate that any observed changes in the early post-blink phase meniscus parameters have to be viewed with caution as they are most likely related to the longitudinal movement of the eye [34] rather than to the tear meniscus formation immediately after a blink that would correspond to the phase of tear film build-up [32,39].
Taking into account the combination of interfering factors related to head, eye and eyelid movements, it is concluded that a single acquisition of tear film meniscus in OCT does not provide reliable estimates of tear meniscus parameters. The clinical utility of dynamic OCT based tear meniscus measurement lies in the necessity of continuous acquisition of post-blink tear film meniscus to provide reliable static estimates of its parameters such as height, depth, and area calculated as averages or medians over a few second post blink interval.