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

Néstor Uribe-Patarroyo; Sahar Hosseinzadeh Kassani; Martin Villiger; Brett E. Bouma

doi:10.1364/OE.26.009081

1. 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–3]. Nonlinear k-sampling and dispersion imbalance are known as two major sources of degradation of the axial resolution in FD-OCT [4–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 path-length 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-in-wavenumber 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 dispersion-compensation 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 k-mapping 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 depth-dependent 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.

2. 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(d). Figure 1(e) illustrates that the full width at half maximum (FWHM) of the peaks increases with depth, showing inconsistent and poor axial resolution without calibration. Extracting the associated fringe signal for a particular reflector position demonstrates the frequency chirp arising from both the nonlinear k-sampling and the dispersion imbalance [Fig. 1(f)].

Fig. 1 Typical configuration of an SS-OCT (a) and SD-OCT (b) system. Image of a mirror at different optical path differences (c), broadening of the peaks due to nonlinearity and dispersion artifacts (d), full width at half maximum (e) of the peaks, and uncorrected, chirped fringe signal (f).

Download Full Size | PDF

In OCT, the spectral interferogram $I_{int}$ of a single sample reflector and the reference arm signal detected by a spectrometer or a photodetector can be expressed as

I_{int} (q, z) = S (q) R_{R} + S (q) R_{s} + 2 S (q) \sqrt{R_{R} R_{S}} \cos [2 f_{k} (q) z + g (q)],

where

k = f_{k} (q)

is wavenumber, and we define

f_{k} (q)

as the wavenumber sampling function.

S (q)

is the spectral envelope of the light source, and

z = z_{R} - z_{S}

is the single-pass OPD between the reference and sample arms.

R_{R}

and

R_{S}

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

{I^{'}}_{i}_{n t} (q, z) = S (q) \sqrt{R_{R} R_{S}} \exp [i φ (q, z)],

where

φ (q, z) = 2 f_{k} (q) z + g (q) .

Equation (3) identifies two contributions to the fringe chirp: nonlinear k-sampling and dispersion imbalance. If we consider two depths

z_{1}

and

z_{2}

, Eq. (3) reads

φ_{1} (q, z_{1}) = 2 f_{k} (q) z_{1} + g (q)

φ_{2} (q, z_{2}) = 2 f_{k} (q) z_{2} + g (q) .

Now, assuming we know

k = f_{k} (q)

we can use its inverse

q = f_{k}^{- 1} (k)

, which we define as the mapping function, to linearize Eqs. (4) and (5)

{\hat{φ}}_{1} (k, z_{1}) = 2 f_{k} [f_{k}^{- 1} (k)] z_{1} + g [f_{k}^{- 1} (k)] = 2 k z_{1} + \hat{g} (k),

{\hat{φ}}_{2} (k, z_{2}) = 2 f_{k} [f_{k}^{- 1} (k)] z_{2} + g [f_{k}^{- 1} (k)] = 2 k z_{2} + \hat{g} (k),

resulting in a common chirp term

\hat{g} (k) = g [f_{k}^{- 1} (k)]

. This means that if we know the correct

f_{k} (q)

, after linearizing the fringes, using

f_{k}^{- 1} (k)

, the residual chirp corresponds to the depth-independent system dispersion

g [f_{k}^{- 1} (k)] = \overset{⌢}{g} (k)

. Throughout this work, we use ^{^} to denote functions resampled to linear k, and ^~ for functions resampled to a general, not necessarily linear, candidate

\tilde{k}

. Now, if we have an incorrect sampling function

\tilde{k} = {\tilde{f}}_{k} (q)

, after remapping using its inverse

{\tilde{f}}_{k}^{- 1} (\tilde{k})

the resulting fringes can be expressed as

{\tilde{φ}}_{1} (\tilde{k}, z_{1}) = 2 f_{k} [{\tilde{f}}_{k}^{- 1} (k)] z_{1} + g [{\tilde{f}}_{k}^{- 1} (k)],

{\tilde{φ}}_{2} (k, z_{2}) = 2 f_{k} [{\tilde{f}}_{k}^{- 1} (k)] z_{2} + g [{\tilde{f}}_{k}^{- 1} (k)],

where

f_{k} [{\tilde{f}}_{k}^{- 1} (k)] z \neq k z

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

f_{k} [{\tilde{f}}_{k}^{- 1} (k)] z = α \tilde{k} z + β O ({\tilde{k}}^{n}) z = α \tilde{k} z + \tilde{γ} (\tilde{k}, z),

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 $\tilde{ψ} (\tilde{k}, z)$ as

\tilde{ψ} (\tilde{k}, z) = \tilde{γ} (\tilde{k}, z) + \tilde{g} (\tilde{k}),

which has a depth dependency for a general

\tilde{k}

, and for

k = \tilde{k}

reduces to

\hat{ψ} (k, z) = \hat{ψ} (k) = \hat{g} (k) .

Therefore, the variance of

\hat{ψ} (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

f_{k} (q)

and

\overset{⌢}{g} (k)

, as an optimization problem.

Assuming a discrete q-sampling with N_s number of samples, we define $f_{k} (q_{j}, \vec{p})$ as an n-order polynomial, and the discrete counterpart to $f_{k} (q)$ , as

f_{k} (q_{j}, \vec{p}) = (N_{s} - 1) γ [c + (\frac{2 q_{j}}{N_{s} - 1} - 1) + \frac{p_{^{2}}}{2} {(\frac{2 q_{j}}{N_{s} - 1} - 1)}^{2} + \frac{p_{3}}{3} {(\frac{2 q_{j}}{N_{s} - 1} - 1)}^{3} + ...],

where q_j is the j-th sample,

\vec{p} = (p_{2}, p_{3}, ...)

is an (n-1)-element vector of polynomial coefficients, and γ and c are defined by boundary conditions. We, therefore, define

f_{k} (q_{j}, \vec{p})

as a relative mapping function, insensitive to the absolute k value in µm⁻¹ of the sampling function. This removes the complexity of determining specific k values and leaves the relation between pixel and µm as an additional calibration step. Note that

f_{k} (q_{j}, \vec{p} = 0)

corresponds to a linear mapping. We consider here the boundary conditions

f_{k} (q = 0, \vec{p}) = 0

and

f_{k} (q = N_{s} - 1, \vec{p}) = N_{s} - 1

. These boundary conditions are arbitrary and depend on how the resampling to linear k is implemented numerically. For our choice, γ and c are given by

γ = \frac{1}{2 + \sum_{l = 2}^{n} \frac{p_{l}}{l} [1 - {(- 1)}^{l}]},

c = 1 - \sum_{l = 2}^{n} {(- 1)}^{l} \frac{p_{l}}{l} .

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 $\tilde{ψ} (\tilde{k}, z)$ is not parametrized. As mentioned above, the variance of $\tilde{ψ} (\tilde{k}, z)$ with respect to z vanishes only for the correct mapping function. For a given $\vec{p}$ , we define this variance as the depth-wise mean square error (MSE) of the residual chirp at sample ${\tilde{k}}_{j}$ as

M S E ({\tilde{k}}_{j}, \vec{p}) = \frac{1}{N_{z} - 1} \sum_{l = 1}^{N_{z}} {(\tilde{ψ} ({\tilde{k}}_{j}, z_{l}, \vec{p}) - \bar{\tilde{ψ}} ({\tilde{k}}_{j}, \vec{p}))}^{2},

where

N_{z}

is number of mirror images at different depths

z_{l}

, and

\bar{\tilde{ψ}}

is the depth average of the residual chirp

\tilde{ψ}

. We define the average MSE across the spectrum as

A M S E (\vec{p}) = \frac{1}{N_{s} - 1} \sum_{j = 1}^{N_{s}} \frac{1}{N_{z} - 1} \sum_{l = 1}^{N_{z}} {(\tilde{ψ} ({\tilde{k}}_{j}, z_{l}, \vec{p}) - \bar{\tilde{ψ}} ({\tilde{k}}_{j}, \vec{p}))}^{2} .

From the analysis above, the $\vec{p}$ that describes the true $f_{k}$ of the system will have a negligible AMSE. Now we pose the search of $\vec{p}$ as an optimization problem: ${\vec{p}}_{\min} = \min_{\vec{p}} {A M S E (\vec{p})}$ , so the AMSE is the objective function to be minimized by searching the space of $\vec{p}$ .

The search process is shown in Fig. 2. The current proposed $\vec{p}$ defines the sampling function $f_{k}$ . This function is tested for invertibility, and if $f_{k}^{- 1}$ exists the process continues. If the sampling function is not invertible, the AMSE is assumed infinite and a new $\vec{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 $f_{k}^{- 1}$ . 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 $\tilde{ψ} (\tilde{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 $\vec{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 $\tilde{ψ} (\tilde{k}, z) ≅ \hat{ψ} (k),$ and therefore the depth-averaged residual chirp $\bar{\tilde{ψ}} (k) = \hat{g} (k)$ provides the dispersion imbalance of the system.

Fig. 2 Flowchart showing optimization algorithm procedure.

Download Full Size | PDF

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 $\vec{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 $\vec{p}$ ; our seed corresponded to $\vec{p}$ = 0 in all cases.

3. 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 $\vec{p} = 0$ . $f_{k}$ 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 $f_{k}$ , and the time required for the optimization procedure. We observed that a 10-order polynomial for $f_{k}$ 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 $f_{k}$ , $d f_{k} / d q$ , has a constant value, as expected for a linear function. Figure 3(b) shows the corresponding mapping function $f_{k}^{- 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², and the initial AMSE was 11.5 rad². 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.

Fig. 3 (a) Initial f_k seed function and its derivative, (b) mapping function and its derivative (c) associated residual chirp at different depths, (d) mean square error of dispersion, (e) optimized f_k and its derivative (f) mapping function for optimized f_k function and its derivative, (g) associated dispersion at different depths after optimization, (h) mean square error of the residual chirp.

Download Full Size | PDF

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⁻⁵, because it provided excellent results for all systems tested. Figure 3(e) shows $f_{k}$ and its derivative, which reveals the small nonlinearity of the polygon-based wavelength-swept laser. Figure 3(f) shows the corresponding mapping function $f_{k}^{- 1}$ that has been used to interpolate the fringe data, and Fig. 3(g) reveals the recovered dispersion $\overset{⌢}{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⁻³ rad²). 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.

Fig. 4 (a) Corrected image of a mirror at different OPDs, (d) PSF at different OPDs, (e) flat FWHM, (f) corrected fringe signal.

Download Full Size | PDF

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 24-order 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 $f_{k}$ function for all four cases, and Fig. 6(a, bottom) the derivative of $f_{k}$ 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 $\vec{p}$ . Figure 6(b) shows the optimum $\vec{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.

Fig. 5 Uncorrected mirror image at different OPDs, corrected mirror image at different OPDs and FWHM respectively for (a)-(c) an FP-based wavelength-swept laser with 25% duty cycle, (d)-(f) an FP-based wavelength-swept laser with 45% duty cycle, (g)-(i) a MEMS scanner-based wavelength-swept laser and (j)-(l) an SD-OCT system.

Download Full Size | PDF

Fig. 6 (a) Nonlinearity of different wavelength-swept lasers: (top) the k sampling function of a FP-based wavelength-swept laser with 25% duty cycle (FP25), a FP-based wavelength-swept laser with 45% duty cycle (FP45), a MEMS scanner-based wavelength-swept laser and a spectral-domain (SD) OCT system. (bottom) The derivative of the k sampling function for the same systems, which more clearly shows the nonlinear behavior. (b) Final $\vec{p}$ for the different systems. Bar height shows the absolute value of each term, light shaded bars correspond to negative terms.

Download Full Size | PDF

4. 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.

5. 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/ .

Funding

National Institutes of Health (NIH) (P41 EB015903); Terumo Corporation.

Acknowledgments

The authors would like to thank Prof. Nima Tabatabaei, York University, Toronto, for providing the SD-OCT calibration data, and SANTEC Corporation for the MEMS scanner-based wavelength-swept laser.

References and links

1. R. Leitgeb, C. Hitzenberger, and A. Fercher, “Performance of Fourier domain vs. time domain optical coherence tomography,” Opt. Express 11(8), 889–894 (2003). [CrossRef] [PubMed]

2. J. F. de Boer, B. Cense, B. H. Park, M. C. Pierce, G. J. Tearney, and B. E. Bouma, “Improved signal-to-noise ratio in spectral-domain compared with time-domain optical coherence tomography,” Opt. Lett. 28(21), 2067–2069 (2003). [CrossRef] [PubMed]

3. M. Choma, M. Sarunic, C. Yang, and J. Izatt, “Sensitivity advantage of swept source and Fourier domain optical coherence tomography,” Opt. Express 11(18), 2183–2189 (2003). [CrossRef] [PubMed]

4. M. Wojtkowski, V. Srinivasan, T. Ko, J. Fujimoto, A. Kowalczyk, and J. Duker, “Ultrahigh-resolution, high-speed, Fourier domain optical coherence tomography and methods for dispersion compensation,” Opt. Express 12(11), 2404–2422 (2004). [CrossRef] [PubMed]

5. S. Yun, G. Tearney, J. de Boer, N. Iftimia, and B. Bouma, “High-speed optical frequency-domain imaging,” Opt. Express 11(22), 2953–2963 (2003). [CrossRef] [PubMed]

6. E. Brinkmeyer and R. Ulrich, “High-resolution OCDR in dispersive waveguide,” Electron. Lett. 26(6), 413–414 (1990). [CrossRef]

7. C. K. Hitzenberger, A. Baumgartner, W. Drexler, and A. F. Fercher, “Dispersion effects in partial coherence interferometry: implications for intraocular ranging,” J. Biomed. Opt. 4(1), 144–151 (1999). [CrossRef] [PubMed]

8. V. M. Gelikonov, G. V. Gelikonov, and P. A. Shilyagin, “Linear-wavenumber spectrometer for high-speed spectral-domain optical coherence tomography,” Opt. Spectrosc. 106(3), 459–465 (2009). [CrossRef]

9. E. Brinkmeyer and U. Glombitza, “High-resolution coherent frequency-domain reflectrometry using continuously tuned laser diodes,” in Optical Fiber Communication Conference, Vol. 4 of 1991 OSA Technical Digest Series (Optical Society of America, 1991), paper WN2.

10. K. Takada, “Fiber-optic frequency encoder for high-resolution OFDR,” IEEE Photonics Technol. Lett. 4(10), 1174–1177 (1992). [CrossRef]

11. S. Yun, G. Tearney, J. de Boer, and B. Bouma, “Removing the depth-degeneracy in optical frequency domain imaging with frequency shifting,” Opt. Express 12(20), 4822–4828 (2004). [CrossRef] [PubMed]

12. C. S. Cheung, M. Spring, and H. Liang, “Ultra-high resolution Fourier domain optical coherence tomography for old master paintings,” Opt. Express 23(8), 10145–10157 (2015). [CrossRef] [PubMed]

13. M. Szkulmowski, S. Tamborski, and M. Wojtkowski, “Spectrometer calibration for spectroscopic Fourier domain optical coherence tomography,” Biomed. Opt. Express 7(12), 5042–5054 (2016). [CrossRef] [PubMed]

14. C. Dorrer, N. Belabas, J.-P. Likforman, and M. Joffre, “Spectral resolution and sampling issues in Fourier-transform spectral interferometry,” J. Opt. Soc. Am. B 17(10), 1795–1802 (2000). [CrossRef]

15. B. Cense, N. Nassif, T. Chen, M. Pierce, S.-H. Yun, B. Park, B. Bouma, G. Tearney, and J. de Boer, “Ultrahigh-resolution high-speed retinal imaging using spectral-domain optical coherence tomography,” Opt. Express 12(11), 2435–2447 (2004). [CrossRef] [PubMed]

16. Y. Yasuno, V. D. Madjarova, S. Makita, M. Akiba, A. Morosawa, C. Chong, T. Sakai, K.-P. Chan, M. Itoh, and T. Yatagai, “Three-dimensional and high-speed swept-source optical coherence tomography for in vivo investigation of human anterior eye segments,” Opt. Express 13(26), 10652–10664 (2005). [CrossRef] [PubMed]

17. S. Vergnole, D. Lévesque, K. Bizheva, and G. Lamouche, “Optimal signal processing of nonlinearity in swept-source and spectral-domain optical coherence tomography,” Appl. Opt. 51(11), 1701–1708 (2012). [CrossRef] [PubMed]

18. R. Huber, M. Wojtkowski, K. Taira, J. Fujimoto, and K. Hsu, “Amplified, frequency swept lasers for frequency domain reflectometry and OCT imaging: design and scaling principles,” Opt. Express 13(9), 3513–3528 (2005). [CrossRef] [PubMed]

19. K. Singh, G. Sharma, and G. J. Tearney, “Estimation and compensation of dispersion for a high-resolution optical coherence tomography system,” J. Opt. 20(2), 025301 (2018). [CrossRef]

20. M. Wojtkowski, T. Bajraszewski, I. Gorczyńska, P. Targowski, A. Kowalczyk, W. Wasilewski, and C. Radzewicz, “Ophthalmic imaging by spectral optical coherence tomography,” Am. J. Ophthalmol. 138(3), 412–419 (2004). [CrossRef] [PubMed]

21. M. J. Todd, “The many facets of linear programming,” Math. Program. 91(3), 417–436 (2002). [CrossRef]

22. T. G. Kolda, R. M. Lewis, and V. Torczon, “Optimization by direct search: new perspectives on some classical and modern methods,” SIAM Rev. 45(3), 385–482 (2003). [CrossRef]

23. C. Jun, M. Villiger, W.-Y. Oh, and B. E. Bouma, “All-fiber wavelength swept ring laser based on Fabry-Perot filter for optical frequency domain imaging,” Opt. Express 22(21), 25805–25814 (2014). [CrossRef] [PubMed]

24. K. Totsuka, K. Isamoto, T. Sakai, A. Morosawa, and C. Chong, “MEMS scanner based swept source laser for optical coherence tomography,” Proc. SPIE 7554, 75542Q1 (2010)

25. R. Huber, M. Wojtkowski, and J. G. Fujimoto, “Fourier domain mode locking (FDML): A new laser operating regime and applications for optical coherence tomography,” Opt. Express 14(8), 3225–3237 (2006). [CrossRef] [PubMed]

26. N. Uribe-Patarroyo, S. H. Kassani, M. Villiger, and B. E. Bouma, “Robust wavenumber and dispersion calibration for Fourier-domain optical coherence tomography,” figshare (2018) [retrieved 15 Jan 2018], https://doi.org/10.6084/m9.figshare.5787840.

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

Abstract

1. Introduction

2. Principle of the method

3. Experimental results

4. Conclusion

5. Supporting material

Funding

Acknowledgments

References and links

Supplementary Material (1)

Cited By

Figures (6)

Equations (16)

Optics Express