Dispersion mapping as a simple postprocessing step for Fourier domain Optical Coherence Tomography data

Optical Coherence Tomography (OCT) was originally conceived as a volumetric imaging method. Quickly, OCT images went beyond structural data and started to provide functional information about an object enabling for example visualization of blood flow or tissue elasticity. Minimal or no need for system alterations make functional OCT techniques useful in performing multimodal imaging, where differently contrasted images are produced in a single examination. We propose a method that further extends the current capabilities of OCT and requires no modifications to the system. Our algorithm provides information about the sample’s Group Velocity Dispersion (GVD) and can be easily applied to any OCT dataset acquired with a Fourier domain system. GVD is calculated from the difference in material’s optical thickness measured from two images obtained for different spectral ranges. Instead of using two separate light sources, we propose to apply a filter-based, numerical procedure that synthesizes two spectra from one broadband spectrum. We discuss the limitations of the method and present GVD values for BK7 and sapphire and ocular media: cornea and aqueous humour of a rat eye. Results corroborate previous measurements using two different light sources.

characterization of highly scattering media. On the other hand, in a method presented in ref. 28 , two light sources are used to generate two A-scans which present a depth profile of a sample at the same lateral point. The two A-scans are then analysed to extract the optical thicknesses of the media under investigation. By comparing these optical thicknesses, GVD can be calculated. This method proved to be successful at determining GVD values in layers as thin as 160 µm and, most recently, allowed dispersion measurements in ocular media such as cornea and aqueous humour in rat eyes 29 .
The experimental configurations in ref. 28 and ref. 29 included two separate, synchronized, swept laser sources at different wavelengths. In these experimental configurations, it is necessary to align light beams coming from the two sources so that they propagate co-linearly towards a sample. Moreover, the fringes are detected by two different photodiodes. These characteristics of the experimental configuration may result in the two images not being perfectly superimposed, and subsequently cause inaccuracies in the measurement.
Here we present a purely post-processing method for GVD-calculation-based material differentiation. This method eliminates alignment issues and can be applied to data from any FdOCT system. We use different glasses of different thicknesses to validate this new approach, and to assess its limitations. Finally, the results for the cornea and aqueous humour of a rat eye are presented.

Materials and Methods
The methods discussed in ref. 28 and ref. 29 and here are conceptually similar to the approach presented in ref. 30 . The authors compare optical thicknesses of layers of water, as well as cornea, aqueous humour and lens in a rat eye, for two different wavelengths. However, in the method of ref. 30 , the difference in optical thickness is used in conjunction with an estimate of the geometrical thickness to determine the change of group refractive index, which then can be easily recalculated to GVD.
In our new approach to the method, two spectra are synthesized by applying two Gaussian filters to the original spectrum. The two filters are centred around two different wavelengths, λ a and λ b (corresponding to two frequencies, ω b and ω a ) in the original spectrum and their total bandwidth doesn't exceed the total spectral bandwidth of the original spectrum (Fig. 1a). The resulting sub-spectra are then independently Fourier transformed to produce two depth profiles at the same lateral position in a layered sample. The optical thickness of a given layer in the sample will be different for these two A-scans, because they correspond to two different spectral regions (Fig. 1b). The walk-off, ∆z ab , is determined from the difference in the optical thicknesses and subsequently used to calculate GVD: where l s is the geometrical thickness of the material, ω ω ω ∆ = − ab b a and corresponds to a wavelength distance a , c -the speed of light in vacuum, and β 2 is the GVD calculated at the midpoint between the two frequencies. Because the midpoint is defined in frequencies, its wavelength counterpart is not an arithmetic average of λ a and λ b , but can be easily calculated with the following formula: λ λ λ λ Since OCT provides information about optical distances, the sample thickness, l s , is unknown. For glass objects used in our study, the precise geometrical thicknesses are given by the manufacturer. However, for the ocular media, the thicknesses have to be estimated since the precise values of n g for the central wavelengths of the filtered sub-spectra are not known. In order to estimate the geometrical thickness of a layer, first, the optical thickness is determined from the A-scan obtained by Fourier transformation of the original spectrum. This optical thickness is then divided by an approximated value of group refractive index: 1.380 for cornea (estimated from ref. 30 ) and 1.341 for aqueous humour (group refractive index of water 31 ).
To validate the method, we use a Spectral OCT system. The light from a supercontinuum light source, Leukos Pegasus, is filtered to have a central wavelength of 875 nm and full width at half maximum (FWHM) of 130 nm. This light beam is sent through a Michelson-Linnik interferometer with achromatic lenses with focal length, f, of 40, 50 and 60 mm in each arm. The signal is then coupled in a single-mode fibre (AFW Technologies FOP830-LB-5-SM800-22), which transmits it to a home-built spectrometer. In the spectrometer, the interfered signal is collimated by a lens (f = 75 mm), dispersed by a grating (1200 lines/mm) and then the resulting spectrum is imaged by the second lens (f = 200 mm) onto the camera (Basler Sprint spL8192-70 km, 8192 pixels, 10 µm pitch). Detected spectra are resampled from wavelengths to wavenumbers and undergo numerical dispersion compensation before they are Fourier transformed 32 .
Limitations of the method. To check the parameters of the system, we used a mirror as an object. The FWHM of the PSF obtained with the original spectrum is 3.5 µm. We filter the original spectrum so that the shorter-wavelength sub-spectrum has a central wavelength of 840 nm and its FWHM is equal to 30 nm, and the longer-wavelength sub-spectrum is centred around 900 nm and its FWHM is also 30 nm. The midpoint is 868.97 nm (see the text under Eq. 1) and the wavelength distance between the peaks, λ ∆ ab − 60 nm (Fig. 1a). The measured PSF's FWHM for the sub-spectra is 13.2 µm and 14.7 µm, respectively. The two PSFs are positioned at the same place in the A-scans for the first interface (Fig. 1b), because the original spectrum is interpolated to convert wavelengths to wavenumbers before filtering and numerically dispersion-compensated before filtering.
The choice of the central wavelengths and the bandwidths that are filtered out from the original spectrum is dictated by several criteria. First, the orders of magnitude of values of β 2 and the thickness of the layer under investigation have to be considered. Too small value of either of these parameters could produce a walk-off that may not be detectable if the spectral separation of the filtered spectra is not large enough (see Eq. 1). The minimum walk-off detectable by an OCT system is related to the resolution of the process of Fourier transformation (FT). The resolution of FT is the minimum difference in frequency which can be distinguished with FT, which for OCT is the minimal distinguishable distance between two peaks in the A-scan. Simple calculations based on basic properties of FT (see Supplementary Information) show that the minimum walk-off, ∆z min , depends on the refractive index of the medium, n, and the total bandwidth of the spectrum, ω ∆ tot , that is Fourier transformed: (2) min tot ω ∆ tot is the total bandwidth of a sub-spectrum. However, in practice, the total spectral bandwidth of the original spectrum should be used as ω ∆ tot in calculations, because the sub-spectra are filtered out from the original spectrum by multiplying it by a window of a Gaussian form, which is non-zero over the entire spectrum and therefore does not limit the sub-spectra's total bandwidth. This information can be used to determine a criterion for the selection of parameters ω ∆ ab and ω ∆ tot for given parameters of a medium: n, β 2 and l s (see Supplementary Information): The accuracy of the technique is further limited by the noise in the system that induces slight fluctuations of the position of the peak in the A-scans over time. These fluctuations produce different optical thickness values and, as a result, slightly different values of β 2 . These values have a Gaussian distribution whose FWHM depends on two parameters that can be controlled in the method: total spectral bandwidth of the sub-spectra and the number of FT points.
To see the influence of these two parameters on the results, we acquired 300 spectra at one point of a 425 µm thick sapphire glass slide. β 2 was calculated for every measured spectrum for FWHM varying from 15 nm to 40 nm for each sub-spectra with a λ ∆ ab of 60 nm and a total number of FT points equal to 524,288. The sample histograms of β 2 values for FWHMs of 20 nm, 30 nm and 40 nm are presented in Fig. 2a-c, respectively. The values of standard deviation of the obtained distributions were plotted as a function of the total bandwidth of each sub-spectra in Fig. 2d.
The error clearly decreases with the bandwidth. The increase in bandwidth improves the axial resolution, which means that the true location of a boundary in an object can be determined more accurately in the A-scan. The situation is equivalent to the energy-time uncertainty principle: the wider the distribution of the energies, i.e. the wider the spectrum of light, the narrower the distribution of time, i.e. the location in the A-scan and, consequently, the determined thicknesses, the walk-off and, eventually, the dispersion.
The minimum detectable walk-off (Eq. 2) determines the minimum thickness of a material for which GVD can be calculated. It decreases as a function of the total spectral bandwidth as presented in Fig. 2e.
Finally, the sub-spectra have to be padded with a sufficient number of zeros to ensure accurate determination of the peaks' maxima after FT. The histograms presenting the distribution of dispersion values for a 425 µm thick sapphire glass for spectra zero-padded to contain 131,072 (2 17 ), 524,288 (2 19 ) and 2,097,152 (2 21 ) points are presented in Fig. 3a-c, respectively.
In the plot in Fig. 3d depicting the error as a function of FT points, we can see that the error increases substantially if the zero-padding is not optimal. Because the number of bins in the histograms in Fig. 3 is fixed, a discrete distribution of β 2 values can be seen as the accuracy in determining the walk-off is one pixel, which corresponds

Results
We chose to perform Fourier transformation on 262,144 (2 18 ) points to ensure a sufficiently accurate determination of the thickness of the layers, which in our case is 0.0364 µm per pixel, for sub-spectra centred around 840 and 900 nm (thus separated by λ ∆ ab = 60 nm) with an FWHM of 30 nm each. The mid-point between the peaks, for which the GVD value is determined, is 868.97 nm (see the text under Eq. 1).  Table 2. Dispersion values and thickness for cornea and aqueous humour of a rat eye for 868.97 nm. β 2 is calculated as a mean of values obtained from 12 A-scans covering lateral distance of 180 µm on the eye. Glass. We started with measuring dispersion values for two double-layered objects: • 425 µm thick sapphire glass on top of 975 µm thick BK7 glass, • 140 µm thick BK7 glass on top of a 975 µm thick BK7 glass. 20 A-scans over a lateral distance of 200 µm on the sample were acquired and analysed in terms of β 2 . Then the dispersion value was calculated as a mean over the dispersion values for each acquired A-scan. The results are presented in Table 1.
The standard deviation of β 2 is larger for thinner samples, because the error of determination of a walk-off is approximately constant -at least when induced only by random fluctuations -and therefore makes the ratio ∆z l / ab s in Equation 1 larger for thinner materials. For 262,144 FT points and assuming these random fluctuations as the only source of error, the standard deviation of the walk-off is around 2 pixels, which is 0.0788 µm and corresponds to 1.80 fs 2 /mm for a 975 µm thick material, 4.14 fs 2 /mm for 425 µm and 12.57 fs 2 /mm for 140 µm.
Rat eye. Next, we went to investigate the dispersion properties of a rat eye measured ex vivo. Figure 4a,b show sample images obtained after Fourier transforming the shorter and longer wavelength regions of the original spectra. The green rectangle indicates a 180 µm by 2.11 mm area of the eye containing 12 A-scans which were used for the dispersion calculation. The most central areas were chosen, so that the error related to refraction on curved surfaces is minimized. This error comes from the fact that two light beams of different wavelengths passing obliquely through a boundary are refracted at different angles and therefore traverse different distances in the sample producing a larger or smaller walk-off.
The geometrical thickness of the cornea and the aqueous humour calculated as a mean over all 12 A-scans is 131.1 ± 4.7 µm and 373.6 ± 4.0 µm, respectively. The expected standard deviation of dispersion value for these layers -assuming only random fluctuations in the system as a source of error -is 13.42 fs 2 /mm for the cornea and 4.72 fs 2 /mm for the aqueous humour, which is why we decided to detect maximum values of peaks that correspond to the boundaries of the two layers (blue dots in Fig. 4c and red dots in Fig. 4d and a sample A-scan in Fig. 4e) in order to determine optical thickness rather than use a segmentation algorithm and then fit a polynomial as this process may introduce additional error. β 2 values for cornea and aqueous humour are presented in Table 2.

Discussion
We have shown that with a simple numerical procedure of filtering the spectra acquired in Fourier domain OCT and comparing optical thicknesses in the resultant A-scans a β 2 value of an imaged object can be determined. Calibration performed on materials with well-defined properties revealed that several parameters have to be optimised for the method to be successful. Firstly, the total spectral bandwidth of the sub-spectra influences the value of minimum measurable walk-off and the size of standard deviation of the results. Next, the number of points used for the FT, if too small, can decrease the accuracy of determining optical thicknesses and substantially increase the error. Finally, the lower limit is set by noise in the system which causes slight fluctuations of the measured optical thicknesses over time and has the most negative impact for thinner materials. The method enabled determination of dispersion values for ocular media: cornea and aqueous humour. The mean value calculated for cornea, 122.7 fs 2 /mm, is approximately 6 times bigger than for water and a little bigger than the value reported in ref. 30 for the similar spectral range, 108.5 fs 2 /mm. The mean dispersion value calculated for aqueous humour, 23.5 fs 2 /mm, is similar to that of water and slightly bigger than the value in ref. 30 for the similar spectral range: 16.8 fs 2 / mm. The same magnitudes of the dispersion value ratios for the cornea and the aqueous humour compared to water were reported in ref. 29 at 1200 nm.
The method of determining β 2 value can be used for monitoring the change of β 2 in layered objects that are semi-transparent to light. Therefore, it excludes highly scattering media where the back surface boundary might not be visible in an A-scan. Nevertheless, this method is more than suitable to detect early signs of diseases that cause a change in dispersion of the tissue such as the eye. However, it must be noted that the layer for which the value of β 2 is monitored must meet the criterion defined in the Equation 3. Moreover, the layer's thickness should be larger than the fluctuations-related error (discussed in the text beneath Eq. 3) of determining the walk-off for that layer.