Dry eye syndrome is a highly prevalent disease of the ocular surface characterized by an instability of the tear film. Traditional methods used for the evaluation of tear film stability are invasive or show limited repeatability. Here we propose a new non-invasive fully automated approach to measure tear film thickness based on spectral domain optical coherence tomography and on an efficient delay estimator. Silicon wafer phantom were used to validate the thickness measurement. The technique was applied in vivo in healthy subjects. Series of tear film thickness maps were generated, allowing for the visualization of tear film dynamics. Our results show that the in vivo central tear film thickness measurements are precise and repeatable with a coefficient of variation of about 0.65% and that repeatable tear film dynamics can be observed. The presented approach could be used in clinical setting to study patients with dry eye disease and monitor their treatments.
© 2015 Optical Society of America
Optical Coherence Tomography (OCT) is a noninvasive imaging modality based on low coherence interferometry . It allows performing real time, high resolution (1 – 10 µm) tomography of ocular tissues and can provide both morphological and functional information. The depth of imaging of OCT is limited by the optical scattering and by the spectral resolution of the detection and is generally about 1 – 3 mm under the surface of the sample. In ophthalmology, OCT is standard of care and is used for the diagnosis and monitoring of several ocular diseases.
Dry eye syndrome (DES) is characterized by a lack of quantity of tears that lubricate the ocular surface . It is caused by decreased tear production and/or increased tear evaporation. DES is one of the most frequent diseases of the ocular surface. Estimates of the prevalence of dry eye range from 7% to 37%, depending on the country, population and definition used . A recent cohort study in the UK showed that around 10% of women have been diagnosed with dry eye and about 20% experienced DES symptoms . DES results in symptoms of tear film instability, ocular discomfort, inflammation of the ocular surface and visual disturbance . It can lead to a damage of the ocular surface and may substantially reduce the quality of life [6,7].
The tear film contains mainly water but also electrolytes, proteins, mucins, sugars, and other water-soluble substances . Its main functions are, on the one hand, the protection, lubrication and nutrition of the ocular surface and, on the other hand, the upkeep of a smooth ocular surface for optimal function of the eye’s optical system. It is composed of 4 layers : the lipid layer (about 40 – 80 nm thick), the aqueous layer (3 – 8 µm thick in healthy subjects), the mucous layer (about 20 – 40 nm thick), and the glycocalyx layer, which consist of meshwork of membrane-bound mucins, which are organic hydrophilic molecules, serving to anchor the aqueous to the ocular surface [9–11]. The lipid layer, although only less than 0.1 μm thick, helps to stabilize the tear film by retarding evaporation .
Several methods have been developed to analyze the tear film and its dynamic in vivo. Traditional methods used for the diagnosis of dry eye include Schirmer’s test [13,14], measurement of break-up time (BUT) [15,16] and ocular surface staining [17,18]. Schirmer’s test comprises the measurement of the wetting of a strip of absorbing paper inserted over the inferior eyelid margin for a fixed time. Ocular surface staining and BUT measurement use fluorescent markers to detect damaged cells, or the apparition of dry spots on the ocular surface, respectively. These traditional methods suffer from considerable limitations . Drawbacks of the Schirmer’s test are its limited repeatability [19,20] and that the test lacks standardization . General limitations of fluorescent marker based methods are their intrinsic invasiveness , subjectivity of the assessment  and induction of reflex tearing [6,21] that modifies the quantity intended to be measured. In addition, ocular staining methods are characterized by poor repeatability and poor sensitivity [19,22].
The measurement of tear film thickness (TFT) is especially tempting in terms of diagnosis and treatment test of DES [23,24], because no matter whether the main cause is related to increased evaporation or decreased tear production, both will finally lead to changes in TFT. Thus, in the past, several attempts have been made to measure TFT. Mishima  estimated the TFT using glass fibers that were put in contact with the ocular surface. Ehlers  developed a method using absorbing paper disks for soaking up human tear film. Prydal et al.  proposed a method based on confocal microscopy. Creech et al.  developed an indirect method based on the measurement of the radius of curvature of the tear meniscus. However, each of the methods has some limitations. In order to diminish surface reflection, confocal microscopy often requires the contact of immersion liquid with the ocular surface, which may cause cornea lesion . Non-contact implementation of confocal microscopy were reported  but are limited by a relatively small field of view and its inability of obtaining depth-resolved tomography. Except the Creech method, none of these methods is suitable for repeated measurements in clinical trials . Apart from the confocal microscopy approach, each method cited in this paragraph is considered to be invasive .
Several noninvasive techniques for the assessment of TFT based on optical interferometry have been developed. The interferences between light reflected at the different layers of the tear film produce interference fringes (bright and dark bands). As analyzed by King-Smith , these interference fringes can be thickness-dependent, wavelength-dependent, or angle-dependent, according to the interferometric method used. Doane  presented a method based on thickness-dependent fringes where the TFT was estimated by counting the fringes in the image plane. One limitation of this method is that for estimating TFT, it is necessary to count the fringes until zero thickness is reached, meaning that one has to wait until the tear film is completely evaporated in order to estimate its thickness [31,32]. While thickness-dependent fringes methods require monochromatic or quasi monochromatic illumination, wavelength-dependent fringes methods require broad bandwidth light sources. Reflectometry is a wavelength-dependent fringes method and is based on the illumination with spatially incoherent thermal light [33–35]. Prydal  developed a method using a spatially coherent illumination (He-Ne laser) that is based on the analysis of angle-dependent fringes. A limitation of the reflectometry approach is its incapability to capture depth cross-sectional visualization of the anterior segment of the eye or visualize directly the tear layers . Spectral domain optical coherence tomography (SD-OCT) can be seen as a particular case of wavelength-dependent fringes method. It differs from reflectometry by three facts: SD-OCT uses spatially coherent illumination allowing a lateral resolution close to diffraction limit, a reference mirror and scanning. By using 2-axis scanning, OCT allows to perform volume imaging of the ocular surface. Several methods based on SD-OCT were reported: TFT was measured indirectly , using ultrahigh-resolution OCT (UHR-OCT) [38–40] and, for contact lenses applications . In addition, a theoretical approach based on decision theory was applied in computer simulations  and on a tear film phantom .
Here we report a new non-contact, noninvasive approach for measuring TFT in vivo in humans at high temporal sampling using a SD-OCT optical setup. The fully automated processing is described. A theoretical analysis of the fundamental limit of the precision of the method is presented. The model and estimation techniques are validated in an in vitro phantom experiment. Precision and accuracy of the in vivo TFT estimation are analyzed. TFT maps measured over time are presented allowing for visualization of the tear film dynamics.
2.1 Experimental setup
The experimental setup used for performing in vivo human ocular front surface tomography is a modified version of an ultrahigh resolution SD-OCT setup reported previously . The setup is composed of a broadband light source, a free space Michelson interferometer and a spectrometer. A schematic diagram is shown in Fig. 1. The light source is a Titanium:Sapphire laser (Integral OCT; Femtolasers Produktions GmbH, Vienna, Austria) providing 170 nm FWHM spectral bandwidth at a central wavelength = 800 nm. The theoretical axial resolution associated with this optical spectrum is 1.2 µm in tissue . The system was designed in order to optimize optical performances for broadband light use. A beam splitter plate (BSN17, Thorlabs GmbH, Dachau/Munich, Germany) with 90/10 splitting ratio was used in the interferometer. The free-space pathway contained two variable neutral density (ND) filters for adjusting the optical power impinging onto the sample and to balance the reference arm signal reaching the spectrometers CCD camera. The collimators used in the interferometer are achromatic fiber collimators (f = 12 mm, Schäfter Kirchhoff GmbH, Hamburg, Germany). For acquiring a volume, we used a two-axis galvanometric scanner (GVS002; Thorlabs GmbH). The beam is focused onto the sample by an OCT scan lens objective (LSM04-BB; Thorlabs GmbH) yielding a transverse resolution of 20 µm on the sample. In the reference arm, a prism pair was used for compensating the group velocity dispersion (GVD) mismatch caused by the optics of this objective.
For human in vivo measurements, the head of the subject is stabilized on a modified slit lamp head rest, and the subject is advised to fixate an internal fixation target. The optical power incident on the subject’s cornea is less than 1.2 mW which is more than 5 times below the maximum permissible exposure (MPE) as specified by the International Electrotechnical Commission (IEC 60825-1)  and by the American National Standard for Safe Use of Lasers (ANSI Z136.1)  to ensure the safety of the optical radiation on the eye.
The spectral signal returning from the interferometer is directed to a custom build spectrometer for detection. The spectrometer is composed of a collimator with focal length f = 100 mm (OZ Optics, Ottawa, Canada) and a holographic transmission diffraction grating with 1200 lines per mm (Wasatch Photonics, Logan, UT, USA). The light emerging from the grating is imaged by means of an objective with a focal length f = 85 mm (ZEISS PLANAR T f/# = 1.4, ZF-IR-I; Carl Zeiss AG, Oberkochen, Germany) onto a high-speed line array CCD camera (e2v EM4CL 2014; Aviva, Essex, UK) with 2048 pixels and maximum readout rate of 70 kHz. The wavelength range of the spectrometer was 680 – 920 nm. The data acquisition and visualization was done in Labview (Labview 2013, National Instruments, Austin, TX, USA), while the post processing was performed in MATLAB (MATLAB ® R2013b, The MathWorks Inc., Natick, MA, USA).
2.2 Model of the tear film OCT signal
The power spectrum at the input of the interferometer and at the output of the interferometer can be linked as following:
where is the magnitude-squared frequency response of the interferometer. The output power spectrum is measured at the spectrometer. Information about the sample’s structure is provided by this filtering of the input power spectrum by .
As the lipid and mucous layers of the tear film have thicknesses smaller than 80 nm, which are far below the theoretical axial resolution of the employed OCT system – which is close to 1.2 µm – these layers cannot be resolved by the OCT imaging system and the tear film is seen as a single layer with two interfaces, as shown in experimental data later in this paper. If the dispersion mismatch between reference and sample arm of the interferometer is perfectly corrected by the adjustable prism pair, the interferometer in Fig. 1 can be reduced to the simplified interferometer layout depicted in Fig. 2, which frequency responsecan be written as:
using the frequency response of the one-dimensional propagation in a weakly dispersing medium Eq. (19), derived in the appendix A.1. Each term of corresponds to a reflector in the interferometer, respectively to the reference mirror, the air-tear interface, and the tear-epithelium interface. Deeper samples reflectors were also observed experimentally. As seen later in this paper, these can be distinguished from the tear film interfaces. Therefore, in sake of clarity, the tear film interfaces are the only samples reflectors considered in this model.
Taking the square modulus of Eq. (2) yields, by neglecting the second order products:
The amplitude factor depends on the splitting ratio s of the beam splitter and the loss factor that is taking the coupling efficiency and the connector losses into account. The electric field reflection coefficient takes, by definition, all losses accumulated during free space propagation (including Fresnel reflection and absorption) into account. The model of the tear film OCT signal presented above is based on several hypothesis. Firstly, in the propagation frequency response [Eq. (19)] used in the model, it is considered that the pulse broadening caused by the group velocity dispersion (GVD) inside the tear film is negligible, an issue that will be discussed later in this paper. Secondly, the second order terms (i.e. with as coefficient) are considered to be negligible in comparison to the first order terms (i.e. with as coefficient). These second order terms are called “autocorrelation artifacts” and correspond to peaks in the coherence function located at the delay differences of each pairs of sample reflectors . In our case, these peaks cannot be observed because their amplitude is below the noise floor. Thirdly, because of the limited axial resolution of the OCT imaging system, the tear film is considered to be seen as a two interfaces structure. It will be shown later in this paper that this simplified model is in good agreement with the in vivo experimental tear film OCT signal data.
Using Eq. (1) and Eq. (3) of the model, it can be seen that the power spectrum at the output of the interferometer is filtered at modulation frequencies corresponding to the delays of the tear film interfaces. Estimating these modulation frequencies in the measured power spectrum allow estimating the delays of each of the two tear film interfaces and can finally be used to determine TFT.
2.3 Tear film thickness
Suppose that the delays and associated with the interfaces 1 and 2 can be estimated [Fig. 2]. The TFT can then be calculated from the delay difference and the group velocity within the tear film by
Using the definition of and Eq. (18), both defined in the appendix A.1, we see that the thickness depends not only on the refractive index at the central wavelength , but also on its first derivative .
2.4 Analogy to frequency estimation problem
Let the vector be
where and is a random vector. The frequency estimation problem consists of estimating the dimensionless frequency from realizations of. Looking back at Eq. (1) and Eq. (3), one recognizes the analogy to the problem of extracting the delays of reflectors from noised OCT power spectrum signal.
2.5 Cramér-Rao lower bound of the delay
Different delay estimators can be applied to the measured OCT signal to estimate a true value, that is a priori unknown. Each estimator will have a specific variance value , which is related to its performance in extracting the delay information under a certain noise condition. The Cramér-Rao lower bound (CRLB) obtains expressions for the minimum possible variance that may be obtained by any unbiased estimator:
This inequality holds for any unbiased estimators of the parameter. The derivation of the CRLB of the delay is shown in appendix A.2. The result writes as
This fundamental bound of the precision of the delay depends only on the OCT axial resolution value and on the signal-to-noise ratio at the peak , as defined in Eq. (29) in the appendix. As shown in appendix A.2, in case of a flat power spectrum embedded into white noise, the constant is close to unity. Other spectral shapes or noise characteristics lead to different values of . If the delay estimator is efficient, it reaches by definition the CRLB. However, other sources of variation of the delay (e.g. motion) have to be eliminated in order to reach this limit.
2.6 Estimator of the delay. Periodogram maximizer
The maximum likelihood (ML) frequency estimation in Gaussian white noise can be shown to be equivalent to the maximization of the periodogram :
The ML delay estimator is then the delay that maximizes of the OCT periodogram over the continuous variable :
The periodogram maximizer delay estimation technique can be shown to be asymptotically efficient, meaning that for large enough N, it reaches the CRLB . An OCT A-scan or axial-scan is a one-dimensional array corresponding to a discretization of the OCT periodogram. One has the following relationship , with where is the number of samples of the A-scan, and . The accuracy of the estimator in Eq. (12), is not limited by the discretization rounding error, but by the SNR and axial resolution only. For performing the continuous periodogram maximization, a standard, derivation-free one-dimensional optimization algorithm is used .
3.1 Validation of the signal model and estimation technique on a phantom
The evaluation of the new approach in vitro was aimed to investigate three major points: firstly, to analyze if the modeled OCT signal of a layered sample as introduced in section 2.2 is in agreement with the measured OCT data of a layered sample; secondly if our method for estimating a delay presented in section 2.6 allows to retrieve accurately a delay, and, finally, if our thickness estimation method is accurate.
For validation of the thickness estimation method, we measured a silicon wafer phantom with layers of known thickness and refractive index with our OCT system and compared the estimated value of the thickness with the reference value.
For fabricating the phantom, a photolithography technique used for semiconductor devices fabrication was used. The phantom was composed of a layer of photoresist (AZ6632 Microchemicals GmbH, Ulm, Germany) placed on a Silicon substrate (Si). The reference thickness of the photoresist layer was obtained with a technique based on white-light interferometry (F20-UVX-thinfilm analyzer, Filmetrics Inc., San Diego, CA, USA) that was considered to give the “true” value.
In Fig. 3(a) a sketch of the fabricated wafer used, composed of a photoresist layer on a silicon (Si) substrate, is shown. The measured power spectrum array was conventionally processed to obtain the OCT A-scan and M-scan . The processing included windowing to diminish the side lobes caused by the light source spectral shape. In Fig. 3(b) a M-scan of the wafer imaged by our UHR-OCT system is presented. An OCT A-scan of the wafer with the ML delay estimates corresponding to optical propagation path of Fig. 3(a) is shown in Fig. 3(c).
We obtained a reference thickness of the photoresist layer (top layer in Fig. 3(a)) of µm using the white light interferometry technique. The delay difference estimated with our method analyzing the OCT signal was , where is the delay of one pixel. The group index of the photoresist at 800 nm is and was calculated using Eq. (18) and the refractive index provided by the manufacturer of the photoresist. The distance in air that corresponds to one pixel was determined via the calibration of the spectrometer using a translation stage and was found to be µm. The two pixel units are linked by , where is the speed of light in vacuum. Using Eq. (4) and these numerical values, we found a thickness of the photoresist layer of µm, in agreement with the reference value . The uncertainty of was calculated using the uncertainty propagation formula. The difference between estimated and reference thickness is about 0.6 nm, or 0.015%, which is within the level of uncertainty. The uncertainty on the estimated value mainly arises from the uncertainty on the distance in air corresponding to one pixel hat was obtained via the translation stage with an uncertainty of 1 nm. This in vitro experiment indicates that our method for estimating a thickness is accurate.
On the OCT M-scan, a third peak (blue curve on Fig. 3(b)) can be observed. We hypothesized that it originates from the double path propagation through the photoresist layer [Fig. 3(a)]. Indeed we found that the delay and amplitude of this peak are compatible with this hypothesis. The discrepancy between the delay difference values and is less than 0.1% if the wafer is precisely aligned perpendicularly to the probe beam of the sample arm. We observed that the discrepancy was increasing up to 15% when the wafer was sufficiently tilted. This can be explained by the refraction of the light inside the photoresist layer, modifying the propagation delay.
3.2 In vivo tear film thickness measurement
Because of eye motion and lower SNR as compared to the in vitro situation, thickness measurement in vivo is more challenging than in vitro phantom measurements. We aimed to investigate the capability of our method to precisely measure the TFT in vivo.
To this end, a series of measurements in healthy subjects was performed and the delay estimator presented in section 2.6 was applied to the measured data. The measurements were carried out in a dimly lit consulting room with a room temperature in the range of 20 – 23°C and humidity in the range of 40 – 60%. The OCT system was operated at an acquisition rate of 70.000 A-scans/s. Each OCT B-scan was composed of 1024 A-scan. One 3-dimensional volume with a size of 4 mm x 4 mm x 1 mm (horizontal x vertical x depth) was acquired within 7 seconds and contained 1024 x 512 x 1024 pixels.
A representative B-scan of the front surface of the cornea measured in vivo is presented in Fig. 4(a). A close up around the apex as depicted in Fig. 4(b) allows to identify the tear film and the corneal epithelium. The depth scale shown in Fig. 4 and Fig. 5 is defined as where is the optical delay, is the speed of light in vacuum, is the tear film group index [Eq. (14)].
As seen in Fig. 4(c), the air-tear and tear-epithelium interfaces are each associated to a peak in the OCT periodogram. Estimating the delay associated with each tear film interface allows determining the TFT, using Eq. (4) and the tear film group velocity. The algorithm used for the delay estimation relies on the maximization of the OCT periodogram and was detailed earlier in this paper. Briefly, the interface delay correspond to the delay of the maximum of the peak. The single A-scan in Fig. 4(c) refers to the location indicated by the white rectangle in Fig. 4(b). The measured TFT value was calculated using Eq. (4), the delay estimator defined in Eq. (12) and the tear film group index indicated in Eq. (14) and equals 5.35 µm. The model of the tear film OCT signal with two interfaces given in Eq. (3) is found in accordance with the experimental tear film OCT data.
In the A-scan depicted in Fig. 4(c), the full-width-at-half-maximum (FWHM) of the air-tear peak is equal to 2.1 µm of tear film equivalent. We see that the measured axial resolution value, defined as the FWHM of the peak, is higher than the expected theoretical axial resolution value of 1.2 µm calculated from the light source’s central wavelength, spectral bandwidth and the OCT axial resolution formula . This is due to the fact that we perform windowing of the power spectrum with Hamming function. This windowing makes it possible to decrease the level of side lobes at the expense of broadening the FWHM of the peaks, i.e. impairing the axial resolution. On the other hand, the OCT axial resolution formula is, strictly speaking, valid only for a Gaussian spectral shape. These facts can explain the difference between the theoretical and measured axial resolution.
The B-scan presented in Figs. 4(a) and 4(b) is composed of A-scan in each column, but are shown in logarithmic scale for better visualization of the whole range of reflectivity. In Fig. 4(b), the periodogram maxima points associated with each tear film interface are shown (blue color for air-tear interface and green color for tear-epithelium interface respectively). In addition, a fourth degree robust polynomial fit of the periodogram maxima points is shown in Fig. 4(b). In order to analyze if the agreement between the fit and the measured points is satisfactory, we calculated two goodness of fit parameters: the coefficient of determination R2 and the root-mean-square error (RMSE). In a typical B-scan, we found R2 = 0.9998 for the air-tear interface (top interface in Fig. 4), and R2 = 0.9992 for the tear-epithelium interface of the tear film. Around the corneal apex (+/− 1 mm), the RMSE was 0.156 for the air-tear interface and 0.353 for the tear-epithelium interface. One explication for higher RMSE of the tear-epithelium interface may be related to roughness of the interface. The non-zero RMSE delay value associated with the air-tear interface is caused by four factors: noise, eye motion relative to the interferometer reference frame, scanning jitters and the geometry of the interface. We compared the measured value of the RMSE with the RMSE values calculated using Eq. (8) from the measured SNR levels around the corneal apex. We observed that the measured RMSE value is more than two orders of magnitude higher than the expected value at the SNR levels and we deduced that the fluctuation of the delay of the interface in a B-scan are almost not caused by noise. The fitting procedure allows finding the average delay of the layer from an ensemble of points and hence is robust against small fluctuations of the delay caused by small movements and scanning jitters.
Concerning the observed SNR levels, as expected, the air-tear interface has a higher SNR than the tear-epithelium interface, which is due to the larger step change in refractive index. The SNR at one pixel of interest was calculated according to Eq. (29). In order to avoid the excess variance due to motion, we calculated the variance of the complex signal where no reflector is present. This does not lead to a biasing of the result as the variance of the noise is uniform . The mean SNR of a layer was calculated in linear scale. The SNR shows a large fluctuation from one point to the other. The mean SNR on the air-tear interface was measured to be 46 dB whereas the maximum SNR was as high as 65 dB. Furthermore, the SNR at the apex or central region was always found to be higher than in the periphery. The maxima of SNR peaks were always observed very close to the apex where the probe beam impinges with an angle close to 90 degrees onto the cornea. The SNR of the tear-epithelium interface, i.e. the cornea front surface, was significantly smaller. Here, the maximum SNR value was measured to be 46 dB, whereas the mean SNR was 32 dB.
Figure 4 shows a measurement with a TFT value around 5.35 µm. However thinner TFT can be observed in some subjects. In Fig. 5, a measurement with a TFT value around 3.1 µm closer to the axial resolution value of the system is shown. For better separation of the two peaks, a different spectral window than the one applied on the measurement shown in Fig. 4 was used. As seen in Fig. 5(c), even if the tear-epithelium interface peak approaches the air-tear interface peak, they can still be distinguished. There is still some margin before reaching the actual resolution limit of the system. This observation indicates that our method can be applied to measure tear films at least as small as 3.1 µm. The horizontal lines in Fig. 4(a) and Fig. 5(a) are due to the background removal operation.
3.3 Central tear film thickness
There can be large variability in the A-scan TFT values, due in particular to the roughness of the tear-epithelium interface. A central TFT value can be associated with each B-scan (see red vertical line segment in Fig. 4(b) and Fig. 5(b)), defined as the central thickness between the polynomial fit associated with each of the tear film interfaces. The post-processing for estimating central TFT is fully automated.
The precision of the central TFT was analyzed by calculating the coefficient of variation, which is defined as the sample standard deviation (SD) divided by the sample mean. The sample was composed of 300 B-scans, each composed of 512 A-scans and acquired in 2.2 s at a fixed position on a healthy subject eye under laboratory conditions. The histogram of the measured value is represented in Fig. 6. The central TFT value was 5.122 ± 0.034 µm (mean ± SD) corresponding to a coefficient of variation of 0.65%. This indicates that the measurement of the central TFT is highly repeatable [Fig. 6].
3.4 Tear film dynamics
Each blink of the eyelids rewets the ocular surface with new tear film. After a new tear film is formed, it undergoes a progressive thinning. The tear film dynamics can be complex and influencing factors include evaporation, gravity, surface tension and tear film breakup effects [49,50].
To investigate if tear film dynamics can be observed using our approach, we acquired 10 volumes just after eye blinking in a healthy subject and calculated TFT maps from each volume. Each volume contained 256 x 256 x 2048 pixels corresponding to 2.5 mm x 2.5 mm x ~1.2 mm and was acquired in 1.5 second, leading to a total acquisition time of 15 seconds. In Figs. 7(b) and 7(c) a series of TFT maps corresponding to a 2.5 mm x 2.5 mm area [Fig. 7(a)] is shown. Each map is numbered time wisely. For the generation of each map, the TFT was calculated in each A-scan using the approach explained earlier in this paper. Afterwards, a median filter of 3x3 kernel size was used to remove salt-and-pepper noise, caused by outlying thicknesses occurring because of bad detection of the tear-epithelium interface. The effect of this filter was analyzed to ensure that it is not causing any bias. The saturation zone on the TFT maps corresponds to the central corneal reflex and is representing a zone where the camera is saturated due to too high optical power caused by an almost perpendicular impact of the beam on the cornea. A possible solution to this limitation is discussed later in this paper. The post-processing for obtaining a series of TFT maps is fully automated and does not require any action of the user.
A progressive thinning of the tear film with time can clearly be observed. After 15 s the total thinning is about 2.5 µm.
To investigate if the observed tear film dynamics are repeatable, a replicate measurement under the same conditions was done directly after the first measurement. The results are shown in Fig. 7(c). When comparing Figs. 7(b) and 7(c), similar tear film dynamics can be observed. This indicates that the measured tear film dynamics after blinking are repeatable in the same subject.
By using a layered phantom of known thickness and group index, we showed that the new approach is accurate for measuring the thickness of a transparent layer in vitro. The measured thickness was found to be in close agreement with the reference thickness. The absolute difference is 0.6 nm, corresponding to a relative difference less than 0.02%.
The central TFT as defined in section 3.3, is a robust measure for tear quantity and can be calculated in each B-scan. The automated algorithm for calculating central TFT is based on the delay estimator defined in Eq. (12) and on robust curve fitting. The precision of the central TFT in vivo was evaluated by calculating the coefficient of variation in a sample of 300 replicate OCT B-scan. A coefficient of variation of 0.65% was found, indicating that the measurement is highly repeatable.
Using ten volumes acquired directly after blinking and analyzing the data according to the above presented approach for TFT estimation, we obtained a series of TFT map allowing to visualize tear film dynamics. To the best of our knowledge, this is the first time that series of two-dimensional TFT maps measured in vivo via OCT are presented. A thinning of the tear film could be observed, which is consistent with the results of King-Smith et al.  and Werkmeister et al. . Furthermore, we showed that repeatable tear film dynamics can be observed. This approach could be used to study dynamic changes in tear film thickness and tear film breakup and has the potential of being an alternative to breakup time measurements (BUT) or non-invasive breakup time (NIBUT) measurements . However, to this end it would be necessary to acquire a larger area on the ocular surface as the image area is limited currently by the low SNR level at larger distance to apex due to the illumination pathway. Recently, Braun et. al.  reported TFT maps using fluorescein and retroillumination imaging of the tear film showing similar thinning pattern as the one presented in the manuscript in hand. King-Smith et al.  presented an interferometric imaging approach based on a modification of Doane’s interferometer  to image the narrowband and broadband interference fringes of the tear film layers over an area of 8x7 mm2 on the ocular surface, allowing to visualize tear film breakup.
The TFT value we measured are consistent with the values obtained by means of other techniques [9,29,38]. Several approaches for measuring TFT based on OCT have already been reported. Chen et al.  reported TFT values of 1.7 ± 1.5 µm and 2.0 ± 1.5 µm. However, the axial resolution in tissue of their system was 3 µm, and TFT values below the axial resolution have to be interpreted with care. Yadav et al.  and Kottaiyan et al.  used an OCT system with axial resolution in tissue of 1.1 µm and presented TFT values of 4.7 ± 1.6 µm in a single B-scan yielding a coefficient of variation of about 38%, details about the algorithm were not reported. Schmoll et al.  used Canny edge detection applied on OCT tomogram with axial resolution in tissue of 1.3 µm and reported central TFT value of 5.1 ± 0.5 µm in a sample of 4 eyes. Werkmeister et al.  used Dijkstra algorithm and reported central TFT of 4.79 ± 0.88 µm in a sample of 26 healthy subjects using a system with axial resolution in tissue of 1.2 µm. In contrast to each of these approaches that used discrete OCT images for the measurement of TFT, our algorithm uses directly the rescaled power spectrum data and relies on a delay estimator which was shown to be efficient , meaning that no other estimator can perform better under a fixed noise condition. Our technique for measuring TFT can be considered as a model-based approach, assuming a model of the measured power spectrum and estimating the parameters of this model. In our model of the OCT signal, the tear film was considered to be one single layer bounded by two interfaces. This simplified model of the tear film was found to be in good agreement with experimental OCT data (cf. Figure 4 (c) Fig. 5(c)). Several other model-based approaches for measuring TFT were reported. King-Smith et al.  and Lu et al.  used fitting of measured reflectance data to estimate TFT, however it did not allow providing en face TFT maps. Huang et al.  presented another model-based approach that was applied in an in vitro phantom.
A limitation of the presented method is that the SNR decreases with increasing distance to the apex. In case the SNR is too low, the tear-epithelium interface cannot be detected. This decrease of the SNR is due to the beam scanning configuration which is almost telecentric. To overcome this limitation, a scanning configuration such that the beam is in every point impinging perpendicularly to the cornea front surface would yield a uniform SNR and thus better precision of the TFT far from the apex. Furthermore, with proper optical power settings, the saturation zone could be avoided. However, the optical design of the scanning system that satisfy this requirement over the whole front surface is complex. On the other hand, the approach is fundamentally limited by the OCT axial resolution. No absolute TFT values smaller than this limit can be measured. This limitation can be overcome using a light source with larger bandwidth yielding smaller OCT axial resolution. This require a proper selection of the wavelength range since a high power spectral density in the visible range can be uncomfortable for subject and disturb the stability of fixation.
4.1 Factors influencing the accuracy of the TFT measurement using OCT
4.1.1 Influence of the probe beam incident angle
If the incident probe beam is not normal to the air-tear interface, a geometrical measurement bias has to be considered. The true TFT is defined as the Hausdorff distance between the incidence point of the probe beam onto the tear film front surface and the front surface of the epithelium. On the other hand, the TFT measured via OCT corresponds to the optical pathway of the probe beam inside the tear film. With the implemented almost telecentric scanning configuration, if the incidence point is distant from the apex of the cornea, the observed TFT appears larger than the true TFT because of a longer propagation distance inside the tear film [Fig. 8(a)]. We modeled this effect using geometrical optics and assuming that the true TFT corresponds to the path along the normal at the incident point on the air-tear interface. Using this assumption, the ratio between the observed TFT and the true TFT can be written asFig. 8(a)]. Using Snell’s law, the value of can be calculated as a function of the distance to apex , the radius of curvature of the cornea, the refractive index of the tear film and the beam angle . The value of is proportional to the voltage of the scanner input signal. In our scanning configuration, the beam angle was measured to be 2.1° at a maximum scanning range of 8 mm x 8 mm in the cornea plane. is the angle between the optical axis of the system and the normal on the cornea at the point of interest.
The result of this model is plotted in Fig. 8(b) for several values of. We see that the TFT value can be biased up to 2.3% at a distance of 2 mm. For this reason, to obtain an accurate TFT measurement, the incidence point has to be close to the cornea apex. With the bias is less than 0.25%. The estimated thickness values can in principle be corrected for this bias if the radius of curvature of the corneal front surface is known. The ideal case would be to have a beam always normal to the corneal front surface to cancel this measurement bias.
4.2 Refractive index and group index of the tear film
4.2.1 Refractive index
Craig et al. measured the refractive index of the tears to be using a digital refractometer . To the best of our knowledge, no measurement of the tear film refractive index over a broad wavelength range is available. The measurement of Craig et al. was done at a single wavelength (589 nm) and under laboratory conditions of pressure and temperature. However, the temperature of the ocular surface ranges between 30 – 35°C and can dynamically vary over time [54–56].
The refractive index of salted water has been investigated in more detail . Empirical formulas with very good accuracy have been developed . These allow for calculating the refractive index of water as a function of wavelength, temperature and salinity. A decrease in refractive index of 0.3% has to be expected when the temperature rises from 20°C to 35°C [Fig. 9]. A change in osmolality can also modify the refractive index, a fact that is relevant because the osmolality of tears can increase after evaporation . Furthermore, an increase in refractive index of about 0.3% can be expected when the salinity changes from 1% to 2% [Fig. 9]. All these parameters have to be taken into account if one aims to measure accurately absolute TFT in vivo.
4.2.2 Group index
It is important to distinguish between group index and refractive index. As can be seen in Eq. (4), for calculating the TFT from the delay difference, the group index at the central wavelength of the light source power spectrum has to be used. No measurement of the tear film group index is available at the central wavelength of our light source. Therefore, we used the group index of water at the osmolality level of tears. As the refractive index of the tear film (blue asterisk in Fig. 9) is not the same as the refractive index of the 1% salted water (blue dashed curve in Fig. 9), we assume an uncertainty of on this value. For the calculation of the TFT, we use:
4.2.3 Influence of the lipid layer
The lipid layer is the outermost layer of the tear film, protecting it from evaporation. The thickness of the lipid layer is estimated to be in the range 40 – 80 nm [33,60,61]. It was measured by two-wavelength techniques or by color look-up-table. The thickness of the lipid layer is far below the OCT axial resolution and can therefore not be resolved. As the refractive index of the lipid is about 1.482  and thus higher than the refractive index of the aqueous, it can lead to a small bias of the estimated TFT values. If one uses only the refractive index of the aqueous to calculate the thickness from the measured delay, the thickness will be slightly overestimated. In case of a lipid layer of 80 nm, which was considered to be the “worst case” when estimating the total TFT, the bias of the TFT is approximately 8 nm.
4.3 Group velocity dispersion
In the model of the OCT signal of the tear film presented in Eq. (3), the group velocity dispersion (GVD) in the tear film was neglected. By definition, the GVD occurs in each media where the second wavelength derivative of the refractive index is non-zero :
The GVD induces a pulse broadening proportional to propagation length. This broadening is symmetrical and does not influence the position of the maximum. Therefore, our delay estimation method based on the maximum detection is robust against GVD in the tear film.
In summary, we developed a novel approach for in vivo TFT measurement. The proposed approach for measuring central TFT is fully automated, has a relative precision of approximately 0.65% and is robust against small eye movements. A new technique for visualization of tear film dynamics was proposed. Repeatable tear film dynamics could be observed. The proposed approach for measuring central TFT is suitable for clinical setting applications and could be used to quantify the dynamic efficacy of DES treatments and for DES diagnosis in combination with other indicators.
A.1 One dimensional propagation of the electric field
In this paragraph, the frequency response modeling the propagation of a collimated beam in a weakly dispersing uniform medium is derived. When propagating through a distance in a weakly dispersing medium, a collimated light pulse travels at the group velocity and the electric field is simply delayed by and multiplied by a phase factor :
where is the group delay, is the group velocity at , and is the group index. in the above equation denotes the phase delay, is the phase velocity at , and is the refractive index. is the central wavelength of the pulse and is the speed of light in vacuum. The group index can be linked to the refractive indexby:
The frequency response corresponding to [Eq. (16)] is:
This frequency response allows linking the complex field Fourier transformation before and after propagation by:
A.2 Cramér-Rao lower bound of the delay
In this section, we derive the CRLB of the delay in a simplified model considering a constant spectral amplitude embedded into white noise. In this model, the power spectrum is represented by the random vector with
where and is a white Gaussian noise random vector of standard deviation . The CRLB of the parameter is [63,64]
where is the signal-to-noise ratio (SNR) associated with . We aimed to obtain the CRLB of the delay. Let us express the same model signal as a function of the delay
where is the optical angular frequency associated with each sample of the rescaled power spectrum, which is given by and, where is the angular frequency separation between adjacent pixels. By identification of Eq. (21) and Eq. (23), one can link the delay and the dimensionless frequency by a scaling
From the domain of and Eq. (24), we obtain the delay range associated with the spectrometer, in accordance with delay range expression found in the literature . In a change of unit of the parameter, the CRLB is multiplied by the square of the scaling factor, as the variance
The total spectral bandwidth of a flat power spectrum is given by
where is the number of pixel. The OCT axial resolution is defined as the FWHM of the peaks in the autocorrelation function . For a spectrum with boxcar envelope, the delay axial resolution is given by 
Let us define the peak SNR as
where is the Discrete Fourier Transform of the random vector, is the modulus of , is the index maximizing the amplitude. It can be shown that the signal-to-noise ratio of and the signal-to-noise ratio at the peak are linked by with the signal model defined by Eq. (21). Using Eq. (27) and Eq. (28), equation (26) can be rewritten as
where is a constant with a value very close to unity . This yields the CRLB of the delay as a function of the delay axial resolution and the SNR at the peak.
Support by the Christian Doppler Laboratory for the Ocular Effects of Thiomers is gratefully acknowledged.
References and links
1. W. Drexler and J. G. Fujimoto, Optical Coherence Tomography: Technology and Applications (Springer, 2008).
2. F. H. Adler, P. L. Kaufman, L. A. Levin, and A. Alm, Adler’s Physiology of the Eye (Elsevier Health Sciences, 2011).
5. M. A. Lemp and G. N. Foulks, “The definition and classification of dry eye disease: report of the Definition and Classification Subcommittee of the International Dry Eye WorkShop (2007),” Ocul. Surf. 5(2), 75–92 (2007). [CrossRef] [PubMed]
6. K. Gumus, C. H. Crockett, K. Rao, E. Yeu, M. P. Weikert, M. Shirayama, S. Hada, and S. C. Pflugfelder, “Noninvasive assessment of tear stability with the tear stability analysis system in tear dysfunction patients,” Invest. Ophthalmol. Vis. Sci. 52(1), 456–461 (2011). [CrossRef] [PubMed]
7. G. Foulks, M. Lemp, J. Jester, J. Sutphin, J. Murube, and G. Novack, “Report of the international dry eye workshop,” Ocul. Surf. 5, 65–203 (2007). [CrossRef]
10. D. M. Albert and F. A. Jakobiec, Principles and Practice of Ophthalmology (Saunders, 2008).
12. J. M. Tiffany, “Physiological functions of the meibomian glands,” Prog. Retin. Eye Res. 14(1), 47–74 (1995). [CrossRef]
13. P. O. Schirmer, “Studien zur Physiologie und Pathologie der Tränenabsonderung und Tränenabfuhr,” Graefes, Arch. Ophthalmol. 56, 197–291 (1903).
14. D. W. Scherz, M. G. Doane, and C. H. Dohlman, “Tear volume in normal eyes and keratoconjunctivitis sicca,” Graefes Arch. Für Klin, Exp. Ophthalmol. 192, 141–150 (1974).
23. S. Kaya, D. Schmidl, L. Schmetterer, K. Napora, A. Unterhuber, V. Aranha dos Santos, C. Baar, G. Garhofer, and R. M. Werkmeister, “Effects of hyaluronic acid on tear film thickness as assesed with high-resolution optical coherence tomography (OCT),” Acta Ophthalmol. (Copenh.). in press.
24. D. Schmidl, L. Schmetterer, K. J. Witkowska, A. Unterhuber, V. Aranha dos Santos, S. Kaya, J. Nepp, C. Baar, P. Rosner, R. M. Werkmeister, and G. Garhofer, “Tear film thickness after treatment with artificial tears in patients with moderate dry eye disease,” Cornea 34(4), 421–426 (2015). [CrossRef] [PubMed]
26. N. Ehlers, “The thickness of the precorneal tear film,” Acta Ophthalmol. (Copenh.) 81, 92–100 (1965).
27. J. I. Prydal and F. W. Campbell, “Study of precorneal tear film thickness and structure by interferometry and confocal microscopy,” Invest. Ophthalmol. Vis. Sci. 33(6), 1996–2005 (1992). [PubMed]
28. J. L. Creech, L. T. Do, I. Fatt, and C. J. Radke, “In vivo tear-film thickness determination and implications for tear-film stability,” Curr. Eye Res. 17(11), 1058–1066 (1998). [CrossRef] [PubMed]
29. T. Schmoll, A. Unterhuber, C. Kolbitsch, T. Le, A. Stingl, and R. Leitgeb, “Precise thickness measurements of Bowman’s layer, epithelium, and tear film,” Optom. Vis. Sci. 89(5), E795–E802 (2012). [CrossRef] [PubMed]
30. R. F. Guthoff, C. Baudouin, and J. Stave, “Principles of confocal in vivo microscopy,” in Atlas of Confocal Laser Scanning In-Vivo Microscopy in Ophthalmology (Springer, 2006), pp. 3–21.
35. P. E. King-Smith, B. A. Fink, N. Fogt, K. K. Nichols, R. M. Hill, and G. S. Wilson, “The thickness of the human precorneal tear film: evidence from reflection spectra,” Invest. Ophthalmol. Vis. Sci. 41(11), 3348–3359 (2000). [PubMed]
36. Q. Chen, J. Wang, A. Tao, M. Shen, S. Jiao, and F. Lu, “Ultrahigh-resolution measurement by optical coherence tomography of dynamic tear film changes on contact lenses,” Invest. Ophthalmol. Vis. Sci. 51(4), 1988–1993 (2010). [CrossRef] [PubMed]
37. J. Wang, D. Fonn, T. L. Simpson, and L. Jones, “Precorneal and pre- and postlens tear film thickness measured indirectly with optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 44(6), 2524–2528 (2003). [CrossRef] [PubMed]
38. R. M. Werkmeister, A. Alex, S. Kaya, A. Unterhuber, B. Hofer, J. Riedl, M. Bronhagl, M. Vietauer, D. Schmidl, T. Schmoll, G. Garhöfer, W. Drexler, R. A. Leitgeb, M. Groeschl, and L. Schmetterer, “Measurement of tear film thickness using ultrahigh-resolution optical coherence tomography,” Invest. Ophthalmol. Vis. Sci. 54(8), 5578–5583 (2013). [CrossRef] [PubMed]
39. R. Yadav, K.-S. Lee, J. P. Rolland, J. M. Zavislan, J. V. Aquavella, and G. Yoon, “Micrometer axial resolution OCT for corneal imaging,” Biomed. Opt. Express 2(11), 3037–3046 (2011). [PubMed]
40. R. Kottaiyan, G. Yoon, Q. Wang, R. Yadav, J. M. Zavislan, and J. V. Aquavella, “Integrated multimodal metrology for objective and noninvasive tear evaluation,” Ocul. Surf. 10(1), 43–50 (2012). [CrossRef] [PubMed]
41. J. Huang, K. S. Lee, E. Clarkson, M. Kupinski, K. L. Maki, D. S. Ross, J. V. Aquavella, and J. P. Rolland, “Phantom study of tear film dynamics with optical coherence tomography and maximum-likelihood estimation,” Opt. Lett. 38(10), 1721–1723 (2013). [CrossRef] [PubMed]
42. J. Huang, Q. Yuan, B. Zhang, K. Xu, P. Tankam, E. Clarkson, M. A. Kupinski, H. B. Hindman, J. V. Aquavella, T. J. Suleski, and J. P. Rolland, “Measurement of a multi-layered tear film phantom using optical coherence tomography and statistical decision theory,” Biomed. Opt. Express 5(12), 4374–4386 (2014). [CrossRef] [PubMed]
43. International Electrotechnical Commission, “Safety of laser products - Part 1: Equipment classification and requirements,” IEC EN 60825–1 Ed. 2, (2001).
44. American National Standards Institute, “American national standard for the safe use of lasers (ANSI Z136. 1),” (2007).
45. B. E. A. Saleh and M. C. Teich, Fundamentals of Photonics, second edition (Wiley-Interscience, 2007).
46. B. G. Quinn and E. J. Hannan, The Estimation and Tracking of Frequency (Cambridge University., 2001).
47. R. P. Brent, Algorithms for Minimization Without Derivatives (Dover Publications, 2013).
50. L. Li, R. J. Braun, K. L. Maki, W. D. Henshaw, and P. E. King-Smith, “Tear film dynamics with evaporation, wetting, and time-dependent flux boundary condition on an eye-shaped domain,” Phys Fluids (1994) 26(5), 052101 (2014). [CrossRef] [PubMed]
51. R. J. Braun, P. E. King-Smith, C. G. Begley, L. Li, and N. R. Gewecke, “Dynamics and function of the tear film in relation to the blink cycle,” Prog. Retin. Eye Res. 45, 132–164 (2015). [CrossRef] [PubMed]
52. P. E. King-Smith, B. A. Fink, J. J. Nichols, K. K. Nichols, and R. M. Hill, “Interferometric imaging of the full thickness of the precorneal tear film,” J. Opt. Soc. Am. A 23(9), 2097–2104 (2006). [CrossRef] [PubMed]
56. M. K. J. Klamann, A.-K. B. Maier, J. Gonnermann, J. P. Klein, and U. Pleyer, “Measurement of dynamic ocular surface temperature in healthy subjects using a new thermography device,” Curr. Eye Res. 37(8), 678–683 (2012). [CrossRef] [PubMed]
57. R. W. Austin and G. Halikas, The Index of Refraction of Seawater (1976).
60. D. R. Korb, The Tear Film: Structure, Function, and Clinical Examination (Elsevier Health Sciences, 2002).
63. S. Kay, “A fast and accurate single frequency estimator,” IEEE Trans. Acoust. Speech Signal Process. 37(12), 1987–1990 (1989). [CrossRef]
64. M. L. Fowler, “Phase-based frequency estimation: a review,” Digit. Signal Process. 12(4), 590–615 (2002). [CrossRef]
65. T. Deleporte, “Quantitative infrared Fourier transform spectroscopy,” Ph.D thesis, Université Libre de Bruxelles (2008).