Modified wavelength scanning interferometry for simultaneous tomography and topography of the cornea with Fourier domain optical coherence tomography

Visual acuity is dependent on corneal shape and size. A minor variation in surface geometry can cause a deformation of corneal geometry, which affects its optical performance. In this work we demonstrate an algorithm for the simultaneous measurement of corneal tomography and topography with a traditional point-scanning Fourier domain optical coherence tomography (FD-OCT) system. A modified wavelength scanning interferometry (mWSI) algorithm enabled topographical evaluation of the surface with nanometer-scale resolution, which is superior to the micrometer-scale resolution of traditional OCT structural imaging. We validated the technique with an optically flat mirror, standard roughness gauges, and atomic force microscopy (AFM). The mirror results show nanometer-scale sensitivity (~3.5 nm), and the mWSI measurements were in good agreement (error ~5%) with the specifications of the roughness comparator and AFM, demonstrating the accuracy of the technique. Following validation, the measurements were made on pig corneas in situ at various artificially controlled intraocular pressures (IOP) and before and after cross-linking (CXL). The results show that the mean surface roughness increased by ~65% after removal of the epithelium in preparation for CXL but did not change as a function of IOP. The demonstrated method could be used for simultaneous measurement of tissue tomography with micrometer-precision and topography with nanometer-precision.


Introduction
Optical coherence tomography (OCT) has emerged as a powerful tool in various fields of biomedical imaging [1]. OCT is based on low coherence interferometry, where the optical path difference (OPD) between the sample and a reference reflector is used to obtain depthresolved images with micrometer-scale resolution. While the major application of OCT has been in ophthalmology, its capability for rapid noninvasive imaging with micrometer-scale resolution has expanded to other fields, such as embryology [2], cardiology [3], dermatology [4], gastroenterology [5], oncology [6], and dentistry [7]. The first realization of OCT was time-domain OCT (TD-OCT) [1], which utilized a scanned reference mirror. The development of Fourier domain OCT (FD-OCT) [8] enabled single-shot depth-wise measurements with superior sensitivity [9] by acquiring spectral information, removing the need for a moving reference mirror. However, the resolution of quantitative structural information provided by OCT is limited to the axial resolution of the system, which is generally on the order of few micrometers.
On the other hand, wavelength scanning interferometry (WSI) can perform topographical surface measurements with much finer resolution than OCT [10,11]. The basic principle of WSI is to shift the phase of the interfered beams by scanning the wavelength instead of altering the optical path difference, such as with traditional phase shifting holography [12]. WSI has the advantage of enhanced stability by overcoming mechanical errors caused by piezoelectric transducers such as hysteresis and nonlinearity, which are inherent to any mechanical scanning method. The key feature of WSI is an acousto-optic tunable filter (AOTF), which selects a specific wavelength from a broadband source, thus producing an interferogram only at that specific wavelength. There have been reports utilizing WSI for topographical evaluation of samples such as integrated circuits and fingerprints via phase reconstruction of wavelength specific interferograms using full-field swept source OCT [13,14]. However, these setups are limited because the AOTF restricts the imaging speed. Additionally, due to the bandwidth limitations of AOTFs, the imaging depth is typically less than a millimeter.
FD-OCT has been extensively used to assess corneal topography [15]. However, only the OCT structural image was used, again restricting the axial displacement sensitivity to a few micrometers. In this paper, we applied a modified WSI (mWSI) algorithm for simultaneous measurement of tomography and topography of the cornea with a traditional point scanning FD-OCT system. Because the demonstrated algorithm does not utilize an AOTF, this method is only restricted in speed by the FD-OCT system and inherits the relatively long imaging range (generally several millimeters). In addition, the mWSI algorithm is applicable for both SS-OCT and SD-OCT systems and equally effective because each A-line is processed independently.
Corneal shape largely determines the optical performance of the eye and even minor variations in surface geometry can alter visual acuity. Irregularities in corneal shape caused by diseases such as keratoconus or scarring due to trauma can severely degrade or impair vision. Most current refractive surgical procedures, such as UV-A/riboflavin corneal collagen crosslinking (CXL), which is used to treat keratoconus [16], work by modifying corneal curvature, and screening for corneal abnormalities prior to these surgeries relies on some form of corneal topography [17][18][19]. To resolve clinically meaningful features, corneal topography measurement techniques must be capable of sub-micrometer axial resolution and micrometer lateral resolution over a 10-15 mm lateral field. Because axial and lateral resolution are decoupled in OCT imaging, it is a promising approach for corneal topography measurements. However, early attempts to use OCT for corneal topography were inferior to existing approaches because imaging resolution (primarily axial) was insufficient. Newer approaches have overcome limitations due to imaging speed and axial resolution has improved [20].
In this work, we describe and evaluate a new approach for OCT-based corneal topography using a point scanning FD-OCT system to achieve nanometer-scale axial topographic resolution. We also present results from topographic measurements of an in situ porcine cornea before and after therapeutic interventions that modify corneal shape.

Experimental procedures
For experimental quantification of the surface topography, we used an optically flat mirror (PF07-03, ThorLabs Inc., NJ, USA), standard surfaces roughness comparators (Models 315 and 316, SPI Microsurf Comparator, MSC Industrial Direct Co., NY, USA). Further validation was performed with an atomic force microscope (AFM) (Bioscope Resolve V, Bruker MA, USA). Please note that AFM validation data are presented in Appendix. Fresh juvenile (4 to 6 months) porcine eyes (Sioux-Preme Packing Co., Sioux City, IA) were used to quantify corneal tomography and topography. The eyes were transported overnight in 4°C saline. Experiments were conducted immediately when the samples were received, which was within 24 hours of enucleation. During the experiments, the whole eye-globes were positioned in a custom holder that was mounted on a stage. After removal of extraneous tissues such as muscles, the eye-globes were cannulated with two needles for artificial intraocular pressure (IOP) control [21]. One needle was connected to a micro-infusion pump, and the other needle was connected to a pressure transducer to form the closed-loop IOP control system. CXL was performed by a technique that mimics the standard "Dresden" protocol [16]. The epithelium was removed with a blunt surgical instrument, and a 0.1% riboflavin solution in 20% dextran was topically applied to the cornea every 5 minutes. After 30 minutes, the cornea was irradiated for 30 minutes (3 mW/cm 2 , 365 nm). During UV irradiation, the riboflavin solution was topically added to the cornea every 5 minutes. The measurements were taken before the treatment, after the epithelium was removed, and immediately after the CXL procedure.

Experimental setup
The SD-OCT system used in these studies has been previously described in detail in [22]. Briefly, the system utilizes a superluminescent diode (SLD) with a central wavelength at 840 nm and −3 dB bandwidth of 50 nm. The source light was split by a 50/50 fiber coupler to a reference arm with a mirror and the sample arm. The light in the sample arm was scanned across the sample surface with a pair of galvanometer mounted mirrors. The light returning from the sample and reference arms interfered, and the interference pattern was captured with a home-built spectrometer. In OCT, the axial resolution (L c ) is determined by two parameters: the central wavelength of the light source (λ 0 ) and the bandwidth (Δλ) by The axial optical resolution of the SD-OCT system was ~6 µm, as measured by the PSF of a mirror, and the transverse resolution was ~8 µm as determined by imaging a resolution target. Each frame took 100 ms with a scanning area of 9.5 mm × 9.5 mm with 2000 × 2000 pixels.

Modified wavelength scanning interferometry algorithm
Considering the lateral dimensions along the surface as the XY plane and depth as the Z direction, the light intensity can be expressed as: Therefore, the dynamic phase change of the interference signal along with the change of wavenumber is given by: The phase difference corresponding to wavenumber shift can be used to obtain the OPD,

( )
, , h x y z , with: There are many methods for phase calculation that can be implemented in WSI. These algorithms include phase demodulation by a lock-in amplifier [22], phase calculation by a seven-point method used in classical phase-shifting interferometry [23], extremum position counting [24], and Fourier transform based techniques [25]. Here, we utilize a Fourier transform based technique because it is fast, accurate, and insensitive to intensity noise [25]. , , , b x y z k define the envelope of the fringe. The envelope slowly varies with respect to k because of the spectral intensity of the light source and spectral response of the detector, which is usually a line-scan CCD camera in an SD-OCT system. The path difference of the interferometer was adjusted to be large enough so that the frequency of the cosine term was greater than the modulation frequency of ( ) , , , a x y z k and ( ) , , , b x y z k for easy separation. Equation (2) can be expressed in exponential form: By taking the Fourier Transform of Eq. (6) with respect to wavenumber, k, we can filter out the first and last terms using a bandpass filter. The first term is the DC signal, and last term is the complex conjugate of the second term, which can both be removed by a bandpass filter. An inverse Fourier transform can then be applied to reconstruct b, which yields the intensity image of the sample under investigation. Accumulating the intensity images frame by frame, we can construct the 3D structural image. Taking the natural logarithm of the newly constructed b value, we can isolate the phase term, where the imaginary term encodes the phase at each pixel: To obtain the phase from the sample surface, the location of the surface was identified from the OCT intensity image. Due to generally high contrast at the air-sample interface at the target surface, a strong discontinuity appears in the intensity of each A-line at the surface. Since, the interference is only from the sample surface in conventional WSI, the phase information directly provides the surface information. However, in FD-OCT, each individual A-line corresponds to a single spectral interference signal. Therefore, a Hanning window, matched to the point spread function of the SD-OCT system, was defined to select the surface information. The Hanning window increased the phase calculation robustness by removing system noise and information not relating to the sample surface. The implementation of Hanning window to extract the raw surface information is an additional modification of the traditional WSI algorithm.
After locating the surface of the sample in each A-line, we performed an inverse Fast Fourier Transform (FFT) on each A-line to reconstruct b as described earlier. The phase information of the reconstructed b can be extracted from the imaginary part of Eq. (7). Since the mathematically calculated phase is limited from -π to π, the phase distribution has discontinuities along wavenumber. After unwrapping the raw phase information extracted from Eq. (7) along wavenumber, the phase follows an inclined straight line. The slope of the line signifies the dynamic change of phase with respect to the change in wavenumber as expressed in Eq. (5). Therefore, the OPD at each scanned position of the mirror surface can be calculated from the slope of the line.
The advantage of this algorithm is that the sectioning depth can be obtained over the full OCT imaging range, which is particularly useful for biological samples that have irregular geometry. Additionally, the slope of the unwrapped phase can be positive or negative, so there is no sign ambiguity in the OPD between the sample and reference arms.

Results and discussion
To test the sensitivity of our algorithm, we first applied it on the interference pattern obtained from a coverslip with the system in a common-path setup (blocking the reference arm of the SD-OCT system). The resulting sensitivity, as quantified by the standard deviation of the OPD of the self-interference, was 3.5 nm at an SNR of 40 dB. We then imaged a simple flat mirror (PF07-03, ThorLabs Inc., NJ, USA). From the intensity data collected from the SD-OCT system, we first tracked the surface of the mirror. Surface tracking was performed on each A-line independently based on the medianfiltered intensity. The intensity plot of a single frame with the tracked surface (in red) is shown in Fig. 1. After tracking the surface of the flat mirror for each frame, we followed the mWSI algorithm to get the OPD for each A-line. The OPDs from each A-scan in a 3D image were then combined.
One limitation of most OCT systems is the presence of the so-called fan (or field) distortion for a non-telecentric scanning system [26]. In general, the effect is related to the scanning architecture of most point-scanning 3-D OCT systems, where perfectly flat surfaces appear curved [27][28][29]. This phenomenon is due to at least two possible effects: 1) the architecture of the scanning system, which is primarily affected by the spatial separation of the mirrors, and 2) design, position, and alignment of the collimating lens in relation to the scanning mirrors. As the result, the flat surface of the mirror appears as a curved surface. To remove the scanning distortion, we subtracted a quadratic curve that was fitted along the curved surface for each B-scan.
A 3D surface plot of the flat mirror without correction is shown in Fig. 2(a). The roughness of the mirror is superimposed on the scanning curvature. As mentioned earlier, we used a quadratic fit to replicate the scanning curvature for each B-scan and subtracted it from the corresponding B-scan. After subtracting the scanning curvature from the 3D surface, the remaining data was used to calculate surface roughness. The 3D examples of surface roughness over a small ~100 µm by ~100 µm region of the mirror at three arbitrary locations are illustrated in Figs. 2(b-d). The statistics of the mirror roughness are provided in Table 1. We discarded values at the periphery of the image due to the low structural image SNR. The effect of SNR on the mWSI algorithm is shown in the Appendix in Fig. 7, where the stability decreases drastically when the OCT image SNR drops below 40 dB. An OCT B-Scan of the optically flat mirror and the corresponding SNR and roughness as calculated by mWSI algorithm are shown in the For further verification of the accuracy of mWSI algorithm, we compared the roughness of known standard surfaces from a roughness comparator (Models 315 and 336, SPI Microsurf Comparator, MSC Industrial Direct Co., NY, USA) to the roughness as calculated by the mWSI algorithm (both measurements in accordance with ISO 4287). We calculated mean roughness value (R a ), which is the average of the absolute values of the profile deviations from the mean line of the roughness. We also calculated the mean roughness depth (R z ), which is the average value of greatest height (peak to valley) of the roughness profiles. The comparison is shown in Table 2. Sample 3D surface roughness plots of three surface comparator samples (N2, N3, and N4 of model 315) are shown in Figs. 3(a-c). The roughness of all comparator samples is plotted in the Appendix in Figs. 9 and 10. The surface comparator results show only a small error in the roughness averaged from all samples (5.1 ± 3.5% for R a and 1.4 ± 11.8% for R z ). Further analysis of the disagreement between mWSI and the specifications of the surface comparator is presented in the Appendix in Fig. 11 with normalized Bland-Altman plots [30]. The average disagreement was 4.4% and 4.5% for R a and R z , respectively. The upper limit of agreement (LoA) for R a was 9.5%, and the lower LoA for R a was −0.8%. Similarly, the upper LoA for R z was 9.5%, and the lower LoA for R z was −0.5%. These results demonstrate that mWSI has a slight positive bias as compared to the manufacturer specifications, but the normalized disagreement decreased as a function of feature size for R a and R z . Imaging of sample N8 of Model 315 also demonstrates the limitation of the technique for rougher samples. If the feature size is above the pixel resolution of the OCT system, the mWSI data is then clipped. Hence, other techniques would need to be utilized in this case. Additional validation was performed with atomic force microscopy (AFM), and the methods and results are provided in the Appendix.  In comparison to existing techniques such as AFM, the presented mWSI method can scan a larger sample area with no contact and much more quickly. However, AFM can provide much finer spatial resolution (angstrom-scale) in contrast to the micrometer (transverse) and nanometer (axial) scale resolution demonstrated by the presented technique. Different groups have proposed surface measurements using optical interferometry, such as absolute distance measurements with a measurement error of 13 nm [31] and large stepped structure with an uncertainty of 45 nm [32]. Moreover, the imaging speed of our method can be further improved by using ultra-fast OCT systems [33].
Next, we performed measurements on an in situ porcine cornea at various conditions. We artificially controlled the eye-globe intraocular pressure (IOP) with a closed-loop feedback system [21]. Imaging was performed at IOPs from 10 to 30 mmHg with 5 mmHg increments under three different treatment conditions: 1) a baseline with no treatment, 2) after removing the epithelium with a blunt surgical instrument, and 3) after traditional CXL [16]. The curvature of the untreated porcine cornea at IOP of 10 mmHg is shown in Fig. 4(a). The surface roughness over a ~100 µm by ~100 µm region from the apex of the porcine cornea at an IOP of 10 mmHg with an intact epithelium is plotted in Fig. 4(b). The roughness of the sample after removing the epithelium and after CXL, also at 10 mmHg IOP, are presented in Figs. 4(c) and (d), respectively, with identical axes and colormap to Fig. 4(b) for better comparison to the sample roughness with the epithelium intact. The roughness of the small areas shown was representative of the total scanned area shown in the figures. Clearly, the roughness as measured by mWSI increased after removal of the epithelium and then further increased in roughness after CXL.
As the presented technique has the advantages of providing both topography and tomography, OCT images of the corneas after each treatment step are also presented. The two-dimensional images of the porcine cornea before the treatment, after removal of epithelium, and after CXL are shown in Figs. 5 (a), (b), and (c), respectively. All the presented images were at an artificially controlled IOP of 10 mmHg. Because of the topical addition of dextran, the cornea thins after CXL [34].  A box plot of the absolute surface roughness amplitude of the porcine cornea at various IOPs and treatment conditions is shown in Fig. 6. The horizontal line is the median, whiskers are the 5th and 95th percentiles, the box is the interquartile range, and the small inscribed box is the mean. The detailed statistical analysis of the surface roughness is given in Table 4 in the Appendix. In Fig. 6, the surface roughness of the cornea increased by ~65% after removal of the epithelium, which is reasonable because of the smooth and homogeneous nature of the corneal epithelium [35]. In addition, the Bowman's membrane is absent in the pig eye [36,37], so removing the epithelium exposed the rougher anterior stroma. Even though the collagen fibrils and proteoglycans are uniformly distributed in deeper layers of the corneal stroma, the collagen fibrils are randomly distributed in the anterior section of stroma [38], which explains the increase in roughness as measured by mWSI after removal of the epithelium. AFM of the human cornea revealed that the diameter of corneal collagen fibrils is 81.73 ± 19.88 nm with mean axial D-periodicity of 68.50 ± 10.3 nm [39]. Since the spatial resolution of collagen fibri collagen fibril Previous w caused the fi roughness of previously, th collagen fibril can also contr axial resoluti corneal surfac As with a pattern. Howe due to variou [44]. Extensiv speckle patte various param applicable to Here, we per (formed by sp correlation be the SD-OCT ils, we are m ls rather than in 6. Box plot of the lobe configuration tions: (untreated) cal instrument, an ng. The whiskers a entral inscribed box ally, we postula e primary caus OP. Even thou surface fibril u our future wo corneal results corneal surface 0]. work has show brils to becom f the corneal his may be du ls, so any micr ribute to the bu on would ena ce only [43]. ll coherence im ever, the speck us parameters, ve research has ern [45][46][47][48]. H meters, such a materials such formed further peckles) and t etween the sur system is a fe measuring the b ndividual fibri absolute surface n at various artifi intact cornea, (ep nd (CXL) after t are the 5th and 95 x is the mean.
ate that the ran se for no mea ugh higher IOP ultrastructure at ork, we will a s as well as hi e topography b wn that CXL in me wavy [42]. he cornea dix). The as well as the corresponding spatial frequencies, was incredibly weak (R 2 <0.05), demonstrating that there was no direct correlation between the speckle in the cornea and the mWSI measurements. However, one has to be cautious by correlating phase measurements using coherence-domain techniques with sample roughness directly as speckle-noise might obscure true surface profile in some situations.
Nevertheless, the proposed algorithm was reasonably accurate as demonstrated by the small error from the manufacturer specifications of a roughness comparator. We also demonstrated surface roughness measurement of an in situ porcine cornea, suggesting the feasibility of applying high resolution topographical imaging of biological samples. Our future work will entail analyzing pathological corneal samples because diseases such as keratoconus alter the corneal fibril arrangement [49][50][51][52]. Moreover, we will measure the surface roughness of primate samples with mWSI to investigate the influence of Bowman's membrane on the surface roughness of the cornea.

Conclusion
In this work, we presented a technique for simultaneous tomography and topography with a modified wavelength scanning interferometry technique that can be utilized with traditional point-scanning Fourier domain optical coherence tomography systems. This technique does not require external synchronization as in a traditional WSI system, thus removing restrictions on imaging speed while also preserving the long imaging depth of OCT systems. We validated the mWSI technique by imaging a mirror and a standard roughness comparator. Our results showed good agreement between the manufacturer provided and mWSI measured roughness, demonstrating the accuracy of the technique. We also imaged an in situ porcine cornea at various conditions demonstrating the applicability of mWSI for characterizing the roughness of tissues.     11. Normalized Bland-Altman plots of (a) R a and (b) R z of the surface roughness comparators, which compares the disagreement between the mWSI-measured roughness and roughness provided by the manufacturer. The blue lines are the lower and upper limits of agreement, and the red lines are the average disagreement between mWSI-measured and manufacturer provided roughness. The gray dashed line is the disagreement as a function of the mean roughness (mWSI and manufacturer specification), and the gray shaded region is the 95% confidence interval for the slope of the disagreement. The correlation and statistical significance of the disagreement as a function of the feature size is indicated in each figure for the corresponding roughness metric. The slope of the normalized disagreement was very small for both R a and R z (−0.003%/nm and −7.7x10 −4 %/nm, respectively).
Speckle is an inherent property of any coherent imaging technique such as OCT. Imaging of surfaces such as glass and metal do not have any obvious speckle pattern in OCT images. Thus, we imaged a sample capable of producing a clear speckle pattern with AFM (Nanoscope Resolve V, Bruker) to further validate mWSI method. A gelatin matrix was deposited on a microscope glass slide. A hole was punched in the gelatin matrix as a landmark for side by side comparison of mWSI and AFM measurements. Figure 12 shows the region imaged by OCT and AFM. Figure 13 shows the comparison between the AFMmeasured and mWSI-measured height profiles and roughness of a region across the glass and gelatin. The results show very good agreement and show the limitations in spatial resolution of the OCT scanning system, which was also unable to resolve the individual collagen fibrils in the cornea samples. A summary of the results is provided in Table 3. The relatively large error in the glass sample may be due to interactions between the AFM tip and the glass.    Here, only the data from the measurements at 10 mmHg IOP are shown. The correlation between the roughness and intensity is indicated in the plot. (b) Spatial frequency power spectra obtained by FFT of the OCT image intensity and mWSI-measured height plotted in (a).
(c) The OCT image intensity power spectra versus mWSI-measured height power spectra for the three samples in (a) and (b). The correlation between the two parameters for each condition is indicated in the figure.

Funding
National Institutes of Health (R01EY022362).