Robust wavenumber and dispersion calibration for Fourier-domain optical coherence tomography

: Many Fourier-domain optical coherence tomography (FD-OCT) systems sample the interference fringes with a non-uniform wavenumber ( k ) interval, introducing a chirp to the signal that depends on the path length difference underlying each fringe. A dispersion imbalance between sample and reference arms also generates a chirp in the fringe signal which, in contrast, is independent of depth. Fringe interpolation to obtain a signal linear in k and compensate dispersion imbalance is critical to achieving bandwidth-limited axial resolution. In this work, we propose an optimization-based algorithm to perform robust and automated calibration of FD-OCT systems, recovering both the interpolation function and the dispersion imbalance. Our technique relies on the fact that the unique function that correctly linearizes the fringe data in k space produces a depth-independent chirp. The calibration procedure requires experimental data corresponding to a single reflector at various depth locations, which can easily be obtained by acquiring data while moving a sample mirror in depth. We have tested both spectral domain OCT and swept source OCT systems with various nonlinearities. Results indicate that the proposed calibration method has excellent performance on a wide range of data sets and enables nearly constant resolution at all imaging depths. An implementation of the algorithm is available online.


Introduction
Optical coherence tomography (OCT) is a non-invasive imaging modality that generates high-resolution, cross-sectional images of the subsurface microstructure of biological tissue. The development of Fourier-domain OCT (FD-OCT) has led to substantial improvements in imaging speed and detection sensitivity, and implementations include spectrometer-based spectral-domain (SD-OCT) and swept-source OCT (SS-OCT), also known as optical frequency domain imaging (OFDI) [1][2][3]. Nonlinear k-sampling and dispersion imbalance are known as two major sources of degradation of the axial resolution in FD-OCT [4][5][6]. Nonlinear k-sampling occurs because the spectrometer angular dispersion, in SD-OCT, and the temporal sampling points, in SS-OCT, are not necessarily evenly spaced in the k-domain. As a result, the detected fringe signals exhibit a chirp that varies as a function of the pathlength difference and, if left uncorrected, produces a significant broadening of the axial point spread function (PSF). In comparison, dispersion imbalance between the sample and the reference arm produces a constant chirp, independent of the sample depth, causing further degradation of the axial resolution [7].
There are two possible approaches to address the nonlinear k-sampling problem: using a special hardware configuration that samples the fringe signal directly in linear k, or acquiring the signal nonlinearly in k and resampling it during post processing to a linear k representation by means of a k-mapping function. For the first solution, a linear-inwavenumber spectrometer using a custom fabricated prism has been proposed for SD-OCT [8]. In SS-OCT, a linear k-sampling method has been developed by incorporating an external clock generator known as k-clock [9,10], which requires advanced acquisition electronics that can handle irregular sampling intervals. Furthermore, the k-clock is not readily compatible with OFDI systems that employ frequency shifting to remove depth-degeneracy [11]. The second solution involves determining the k-mapping function that describes the actual k value being sampled as a function of camera pixel (SD-OCT) or time (SS-OCT), either on the fly or in a separate calibration procedure [11]. Traditional k-mapping calibration of the FD-OCT interferometric signal is performed by imaging a mirror at one or two different axial positions and optimizing for either minimum width of the peak or maximum peak signal intensity. This procedure does not guarantee optimum performance at all imaging depths, because optimizing such metrics can lead to artificial mapping and dispersion curves which compensate each other only at specific depth positions. In SD-OCT, calibration of the spectrometer is commonly performed using standard narrow spectral line emission from a calibration lamp, which requires additional hardware. The pixel position versus wavelength data is then fitted with the grating equation or a polynomial [12]. Recently, a new method to calibrate a spectral OCT spectrometer using the Doppler frequency of a sample moving at a constant speed has been introduced [13].
Once the k-mapping function is known, the process of resampling the interferometric fringes to a linear k representation is generally done numerically in post processing. The most common approach is remapping the spectral data to equal sampling intervals in k-space by interpolation prior to discrete Fourier transform (DFT) [12]. In addition to the commonly used linear interpolation [14], cubic spline interpolation has also been used but is computationally intensive [15]. Interpolation via DFT zero-padding is a computationally more efficient strategy to improve the accuracy of the resampling [11]. Resampling through convolution with a Kaiser-Bessel window followed by DFT has also been demonstrated to suppress interpolation artifacts at competitive computational cost [16]. A nearest-neighbor-check algorithm was introduced to minimize processing time at the cost of lower resampling accuracy [17].
With respect to dispersion imbalance, perfect hardware-based dispersion compensation is challenging to achieve, and therefore compensation in post processing is frequently performed. It has been shown that bandwidth-limited resolution can be recovered by numerically compensating the dispersion, therefore reducing the need to use dispersioncompensation blocks as a way to perform hardware-based dispersion balance [18]. The exact dispersion mismatch is typically unknown and needs to be calibrated, which can be easily accomplished by finding the chirp in the fringe of a mirror reflection, once this fringe is known as a function of linear k [19]. FD-OCT is particularly well suited for numerical dispersion compensation as it provides direct access to the spectral interferogram. An alternative approach for numerical dispersion compensation is to model the dispersion with a polynomial and optimize an image sharpness metric by varying the polynomial coefficients [4,20]. However, optimizing based on a metric related to PSF width or sharpness has only been used successfully in systems close to calibration or without large nonlinearities in k sampling. Our own experience indicates that when the k-sampling is unknown and highly nonlinear, PSF metrics do not evolve through a monotonic path to the global minimum during optimization, because the PSF shape goes through many shape transformations with varying sidelobe magnitude. Therefore, this kind of optimization leads to local minima.
Here, we present a novel optimization-based calibration method that provides both the kmapping function and the dispersion imbalance of the system. The required calibration data set is obtained by recording the signal of a single reflector positioned at various axial locations. As stated above, the dispersion imbalance of the system corresponds to a residual chirp in the properly mapped linear-in-k fringes. The key insight behind our algorithm is the fact that this dispersion chirp is identical for the fringes of all depths only if they have been properly resampled in linear k. An incorrect k resampling will yield an erroneous depthdependent dispersion chirp. We note here that our goal is to compensate for system dispersion imbalance, not sample-induced dispersion that indeed would lead to a depth-dependent dispersion chirp. Our method does not require any prior knowledge of the mapping or dispersion of the system, and the calibration data set can be easily acquired by recording the signal of a sample reflector while sweeping the path length difference. By focusing on the variation across depth of the residual chirp after resampling by a candidate k-mapping function, we avoid the drawbacks of PSF-based optimization procedures. We also focus on a very general calibration technique that can be used in any OCT system with no prior knowledge, with robust performance and without the need for fine tuning parameters for each individual system. We have tested different types of wavelength-swept lasers, ranging from the almost-linear wavenumber tuning of a polygon-based laser to very nonlinear sweeps of Fabry-Perot (FP)-based and MEMS scanner-based lasers. In addition, we have also applied this technique to calibrate an SD-OCT system. The results show that the proposed calibration method has an excellent and consistent performance among different configurations with varying degrees of nonlinearities and dispersion imbalance.

Principle of the method
We acquired calibration data sets in which a mirror was placed in the sample arm and fringes were acquired while scanning across different axial positions. We moved the reference arm mirror instead of the sample mirror for simplicity, to avoid signal attenuation due to the shallower depth of field of the focal spot in the sample arm. The signal of the sample mirror was attenuated in order to avoid saturation artifacts. Figure 1(a) shows the schematic of the OFDI system we used to test the wavelength-swept lasers, and Fig. 1(b) shows the schematic of the SD-OCT system. Figure 1(c) presents the tomogram of such a calibration data set showing the mirror at different optical path differences (OPDs). The corresponding PSF as a function of depth is displayed in Fig. 1  In OCT, the spectral interferogram int I of a single sample reflector and the reference arm signal detected by a spectrometer or a photodetector can be expressed as is wavenumber, and we define ( ) k f q as the wavenumber sampling function. ( ) S q is the spectral envelope of the light source, and R S z z z = − is the single-pass OPD between the reference and sample arms. R R and S R are the reference and sample arm reflectivities, respectively. ( ) g q is the dispersion mismatch between the two optical paths. Here, q is time-sample number for SS-OCT and pixel number for SD-OCT. The first two terms of Eq. (1) can be removed or ignored by proper suppression of the DC component of the signal. The desired information is thus encoded in the last interference term. To obtain the analytic complex fringe signal we use the Hilbert transform, taking the DFT of the original fringe signal and suppressing one of the two mirror terms by setting all negative (or positive) frequencies to zero, followed by transforming back to the original domain by inverse DFT. The resultant complex fringe signal is then given by Equation (3) identifies two contributions to the fringe chirp: nonlinear k-sampling and dispersion imbalance. If we consider two depths 1 z and 2 z , Eq. (3) reads , which we define as the mapping function, to linearize Eqs. (4) and (5) resulting in a common chirp term This means that if we know the correct ( ) . Throughout this work, we use ^ to denote functions resampled to linear k, and ~ for functions resampled to a general, not necessarily linear, candidate k  . Now, if we have an incorrect sampling function ( ) after remapping using its inverse 1 ( ) k f k −   the resulting fringes can be expressed as and is therefore, in general, a function not linear in k. This function can be viewed as the sum of a linear term and a non-linear term in k where n ≥ 2. Clearly, the non-linear term produces a different chirp at each sample depth, and is easily identified as an erroneous depth-dependent dispersion. We now define the non-linear-in-k component of ( , which has a depth dependency for a general k  , and for k k =  reduces to ˆˆ( , ) ( ) ( ). k z k g k ψ ψ = = Therefore, the variance of ˆ( , ) k z ψ with respect to z vanishes for the correct mapping function. This property of the calibration data set makes it possible to express the calibration of the system, that is finding ( ) k f q and ( ) g k  , as an optimization problem.
Assuming a discrete q-sampling with N s number of samples, we define ( , ) where q j is the j-th sample, 2 3 ( , ,...) is an (n-1)-element vector of polynomial coefficients, and γ and c are defined by boundary conditions. We, therefore, define ( , ) It is important to note that the polynomial parametrization only applies to the k sampling function and is not related to the dispersion of the system. The dispersion ( , ) k z ψ   is not parametrized. As mentioned above, the variance of ( , ) k z ψ   with respect to z vanishes only for the correct mapping function. For a given p  , we define this variance as the depth-wise mean square error (MSE) of the residual chirp at sample j k  as ( ) where z N is number of mirror images at different depths l z , and ψ is the depth average of the residual chirp ψ . We define the average MSE across the spectrum as ( ) From the analysis above, the p  that describes the true k f of the system will have a negligible AMSE. Now we pose the search of p  as an optimization problem:  , so the AMSE is the objective function to be minimized by searching the space of p  .
The search process is shown in Fig. 2. The current proposed p  defines the sampling function k f . This function is tested for invertibility, and if 1 k f − exists the process continues. If the sampling function is not invertible, the AMSE is assumed infinite and a new p  is proposed. To resample the fringes, we first use interpolation of the complex fringes using DFT zero-padding with an oversampling factor of 8, and then use linear interpolation to map the fringes according to 1 k f − . The next step is to analyze the residual chirp. To this end, the phase of the fringes is extracted, the linear component subtracted, and the nonlinear phase is unwrapped to generate phase curves. To accurately estimate the linear component, the first moment of the DFT of the complex fringes (that is, of the tomogram) is calculated, which allows us to calculate the corresponding linear phase slope and subtract it from the fringe phase. The resulting fringe phase is unwrapped and fitted with a linear polynomial to further subtract any remaining linear component. We found this procedure significantly more robust than directly unwrapping the fringe phase, especially for fringes with frequencies close to the Nyquist limit. At this point, the remaining nonlinear phase corresponds to ( , ) k z ψ   . Next, the AMSE is calculated according to Eq. (13), which represents the variance of the dispersion at all depths. Based on the AMSE, the optimization routine selects a new candidate p  for the next iteration, and subsequent iterations reduce the value of the AMSE. The optimization stops when the change in the AMSE is lower than a defined tolerance. At this point the residual chirp is mostly independent of depth ( , ) and therefore the depth- provides the dispersion imbalance of the system. The algorithm used by the optimization procedure can be a generic one and is a matter of choice, hence, we do not implement a custom routine. The specifics on how p  is updated in response to the AMSE (Fig. 2) depend upon the optimization routine. In our case, we tested two different approaches. First, we applied a local optimizer based on the Nelder-Mead simplex method [21], readily available in Matlab as fminsearch. This simplex-based method demonstrated good performance across all SS-OCT lasers that we tested. However, when the wavenumber sampling was highly nonlinear this method provided suboptimal performance.
The second approach was the global optimizer pattern search [22], available in Matlab as patternsearch, which showed excellent performance for all systems. The simplex algorithm, as with all local optimizers, is sensitive to the seed for p  ; our seed corresponded to p  = 0 in all cases.

Experimental results
We first show in detail our technique for the polygon-based wavelength-swept laser described in Sec. 2. The polygon laser has a 3dB-bandwidth of 103 nm centered at 1290 nm and a 54 kHz repetition rate. Figures 3(a)-3(d) show the first iteration of the optimization, which makes use of the seed 0 p =  . k f and its derivative with respect to q are shown in Fig. 3(a). The order of the polynomial will be a balance between the accuracy of the parametrization of k f , and the time required for the optimization procedure. We observed that a 10-order polynomial for k f offered the best balance for most of the wavelength-swept lasers we tested, and was used for the data shown in Fig. 3. The derivative of k f , k df dq , has a constant value, as expected for a linear function. Figure 3(b) shows the corresponding mapping function 1 Fig. 3(c) indicates the residual chirp attributed to the dispersion of the system at different depths, and Fig. 3(d) displays the MSE of the dispersion as a function of sample number. The MSE function has large values reaching 60 rad 2 , and the initial AMSE was 11.5 rad 2 . To ease visualization, we normalized the sampling and mapping functions from 0 to 1. This allows us to easily compare the calibration of different systems. Figures 3(e)-3(h) shows the final result after optimization using simplex search. Because the AMSE is normalized by the number of samples, its value can be used to compare the calibration quality between different systems. We used a tolerance value of 10 −5 , because it provided excellent results for all systems tested. Figure 3(e) shows k f and its derivative, which reveals the small nonlinearity of the polygon-based wavelength-swept laser. Figure 3(f) shows the corresponding mapping function 1 k f − that has been used to interpolate the fringe data, and Fig. 3(g) reveals the recovered dispersion ( ) g k  . As can be seen, the dispersion has a nearly constant value across depth suggesting that the optimized parameters indeed accurately describe the mapping of the system. Figure 3(h) shows the MSE of the dispersion after optimization, which has a very small average value (AMSE = 2.5 × 10 −3 rad 2 ). Here, the total processing time was 11 min for the simplex method and 20 min with pattern search using a standard PC workstation. However, this time corresponds to the processing time for a system with no prior knowledge of its calibration parameters: when concerned with recalibration of lasers which drift over time, such as FP-based and MEMS scanner-based lasers, processing time can in principle be made shorter by using the last known calibration as starting point in the optimization process. After the optimization, the calibration data set was resampled and dispersion compensated with the retrieved mapping and system dispersion. Figure 4(a) displays the calibration mirror image at different OPDs. Figures 4(b) and 4(c) show the PSFs at different depths and their FWHM, respectively. The FWHM is invariant along depth and the experimental value of FWHM is very close to the calculated transform-limited value. Figure 4(d) depicts the correctly interpolated and dispersion-compensated fringe corresponding to the raw fringe in Fig. 1(f), with no residual chirp. Now we show the calibration results for four different cases, spanning from low, to high non-linear k sampling. For the case of low non-linear k sampling, we used a FP-based wavelength-swept laser in which the FP filter was driven by a sinusoidal signal at the resonant frequency of the filter. The intrinsic nonlinear sweeping arising from the sinusoidal modulation of the FP filter was limited by modulating the semiconductor optical amplifier (SOA) of the ring laser to have a duty cycle of 25% [23], in order to only utilize the most linear part of the sweep. For the case with medium nonlinearity, we increased the duty cycle to 45%. In both cases, the central wavelength of the sweep was 1280 nm. The 3dB-bandwidth was 112 nm for a 25% duty cycle (240 kHz repetition rate) and 150 nm full-bandwidth for a 45% duty cycle at 120 KHz. High non-linear k sampling was obtained with a MEMS scanner technology (Santec Corp.) [24] in which the lasing wavelength was tuned sinusoidally at the resonant frequency of the MEMS mirror. The MEMS-based laser had 105 nm full-bandwidth centered at 1310 nm and 100 KHz repetition rate. In the fourth case, we applied the proposed optimization algorithm to a custom-built spectral domain OCT system with 105 nm 3dB bandwidth and 1315 nm central wavelength. We calibrated the four cases with our optimization-based technique and obtained excellent results in all instances. We show in Figs. 5(a), 5(d), 5(g) and 5(j) the uncalibrated images of the calibration mirror located at different depth positions for all systems. Figures 5(b), 5(e), 5(h) and 5(k) show the calibrated mirror images for all light sources. Figures 5(c), 5(f), 5(i) and 5(l) show that the FWHM of all systems is very nearly constant with depth, and very close to the transform-limited value. A 10-order polynomial for f k and simplex search was used for all SS-OCT systems, except for the FP-based laser with 45% duty cycle: although this laser has a moderately linear sweep at the center of the spectrum, it becomes highly nonlinear toward the edges of the spectrum. A 32-order polynomial and the pattern search method was used for this source. Although this k sampling could be described by a set of a few low and high orders, we would need a priori information on the specific nonlinearity to select these orders, which we were striving to avoid. We found that pattern search was more robust when the polynomial order became very high (>24). The k-sampling behavior of the SD-OCT system is similar to that of the FP-based laser with 45% duty cycle, therefore a 24order of polynomial with simplex search was used. Note that the transform-limited PSF was computed from the shape of each light source spectrum. No window was applied prior to DFT in any case. The FWHM metric does not take into consideration possible sidelobes in the PSF. Therefore, a PSF with a shape close to the transformed-limited PSF, with slightly larger sidelobes, will exhibit a FWHM below the transform-limited FWHM. This explains why, in some cases, the FWHM seen in Fig. (5) is below the transform-limited value. In the case of Fig. 5(f), we found that the system exhibited a significant amount of polarization mode dispersion, which affected the FWHM results away from the center of the imaging range. Figure 6(a, top) shows the normalized optimum k f function for all four cases, and Fig. 6(a, bottom) the derivative of k f to more clearly see the differences in nonlinearity. As it can be seen, the MEMS scanner-based wavelength swept laser exhibited the largest non-linear k sampling among all investigated systems, and the level of nonlinearity was similar throughout the sweep range. In contrast, the FP-based laser with 45% duty cycle and the SD-OCT system exhibited strong nonlinearities only at the edges of the sweep range, which explained the need for higher order terms in p  . Figure 6(b) shows the optimum p  for all systems. It is evident that the FP-based laser with 45% duty cycle and the SD-OCT system have the most weight in a few low-order terms and most high-order terms.

Conclusion
Calibration of non-uniform k mapping and dispersion imbalance in OCT systems is necessary to achieve bandwidth-limited axial resolution at all depths. Despite this, there is a general lack of tools in the literature designed to perform simple and reliable calibrations of general OCT systems. In this work, we presented a calibration method that has the potential to fill this gap. We have tested this technique on different wavelength-swept lasers with varying nonlinearities, as well as on an SD-OCT system. The results show excellent performance of this calibration method in a wide variety of systems with different degrees of non-linear k sampling. In all cases, we achieved an axial resolution close to the bandwidth-limited resolution at all imaged depths. This technique requires a calibration data set that is simple to acquire, and which can be easily automated in, e.g., clinical systems with no user intervention. The robustness of the optimization-based approach simplifies the analysis of the calibration data set, enabling a fully automated calibration of OCT systems. This may be particularly useful for systems that are highly sensitive to temperature variations [25], and can be envisioned as an automated calibration procedure after each system startup and prior to imaging.

Supporting material
A Matlab implementation of our algorithm, together with exemplary data, is in the supplementary materials, Code 1 (Ref [26]), and also available on http://octresearch.org/ .