Constrained polynomial fit-based $k$-domain interpolation in swept-source optical coherence tomography

1. Introduction

Fourier-domain optical coherence tomography (FD-OCT), which is promising for in vivo high-resolution biomedical imaging, has replaced time-domain OCT as the mainstream technology due to its prominent improvement in speed and signal-to-noise ratio (SNR).^1,2,3 In FD-OCT, depth information is acquired by Fourier transform of spectrum as wavenumber-depth ($k$–$z$) is a pair of Fourier transformation. Nonuniform discrete Fourier transform (NDFT) reconstructs depth profiles by dot products between spectrum and a Vander-Monde matrix consisting of nonuniform $k$.^4,5 Although NDFT can obtain precise results, equivalent to that acquired by the discrete Fourier transform, its algorithm is of high complexity and time-consuming. In vivo dynamic monitoring of biological tissue, such as small blood capillary network in the cortex and cerebral fluid,^6,7,8 calls for real-time imaging where massive raw spectrum data are prepared to be processed and analyzed. Graphics processing unit accelerated nonuniform fast Fourier transform is applied to ultrahigh-speed, real-time FD-OCT.⁹ Moreover, Fast Fourier transform (FFT) based on uniform sampling is more widely used in data processing.^10,11 However, in these cases, spectral data are necessary to be linear-in-$k$ in prior; otherwise, the resolution and precision of transformed depth information will be degraded to a large extent.¹²

Two implementations in FD-OCT can acquire linear-in-$k$ signals with unique and elaborate hardware. In Spectral domain optical coherence tomography (SD-OCT), $k$-space spectrometer is designed for linear-in-$k$spectral data acquisition,^13,14,15 while linear-in-$k$ swept laser can ensure uniform sampling in swept-source optical coherence tomography (SS-OCT).^16,17,18 However, these devices have exquisite structures that require a strict and regular calibration. Therefore, resampling algorithm becomes more widely researched and adopted.^19,20,21 In SS-OCT, real-time $k$ calibration refrains the accuracy of resampling for the fact that central wavelength shift exists in each swept period of the swept-source.²² Mach–Zehnder interferometer (MZI) is connected with swept-source to provide spectral phase for real-time $k$ calibration and subsequent time-domain spline interpolation (TDSI).²³ The precision of TDSI is proved to hinge upon order of time versus wavenumber polynomial fit, so $k$-domain spline interpolation is proposed to release processing burden for improved time-efficiency.²⁴

In this paper, a $k$-domain spline interpolation method based on constrained polynomial fit (CPF-KDSI) is proposed. Due to a higher SNR in central region than its surroundings of the MZI signal, we add the constraint condition of passing through the central spectral point to the polynomial fitting of the spectral phase for wavelength calibration. We acquire the point-spread-function (PSF) of the system and the reconstructed OCT images using CPF-KDSI, and compare them to that acquired by conventional KDSI.

2. Method

2.1. Conventional KDSI

In SS-OCT, interference signal is sampled linear-in-time, so an extra $k$ calibration signal in sync is needed for resampling of interference spectrum. An MZI unit can be linked to the swept-source to provide a real-time calibration signal. The MZI calibration signal without regard of noise can be expressed as

where $t_i$ is the timestamp of the $i$th sampling point, $R$ is a scale factor related to transmittances of two arms of the MZI, $S(t_i)$ is the source power spectral intensity, $\phi (t_i)$ is the phase term of the MZI calibration signal, $k(t_i)$ is the wavenumber, and $d$ is the optical path difference (OPD) between two arms of the MZI. Phase term of Eq. (1) can be extracted by Hilbert transform as

One should note that the extracted phase is wrapped in [$-\pi ,\pi$] because of “$2\pi$ ambiguity”.²⁵ Wrapped phase can be unwrapped by adding integral multiples of $2\pi$ at jump points under the fact that the spectral phase difference between adjacent sample points is much smaller than $2\pi$. Conventional KDSI directly linearizes the phase distribution by making use of unwrapped phases of the first and the $N$th sample points

According to Eq. (2), the phase varies linearly with $k$. Therefore, the resampling interference signal using uniform phase coordinates, as shown in Eq. (3), results in the same linear-in-$k$ spectral data with uniform $k$ coordinates.

2.2. KDSI with constrained polynomial fit

Unfortunately, there is relatively a large phase error at the beginning and the end of unwrapped phase in practice. One reason is that the MZI signal is lower in intensity at edges than around central region and a larger phase error occurs in low SNR area where noise exists. Another issue is the introduced error by Hilbert transform at edges of the extracting phases due to the finite length of the sampled signal.²⁶

In CPF-KDSI, the least-squares polynomial fit of phase versus time is constrained by passing through the sample point $N$/2 (($N+1)$/2 when $N$ is odd)

where ${\phi }_0$ is the unwrapped phase after the Hilbert transform. Phase distribution from the fitting curve in Eq. (4) can be linearized by

$${\phi }_F^{\prime }(i)=\frac{{\phi }_F(t_N)-{\phi }_F(t_1)}{N-1}\cdot (i-1)+{\phi }_F(t_1),i=1,2,\dots ,N.$$(5)

Figure 1 illustrates how conventional KDSI and CPF-KDSI obtain uniform phase coordinates for interpolation. Spectral phase from MZI intrinsically varies smoothly without abrupt fluctuations, so the phase from fitting curve is more proximity to the actual phase than the unwrapped one. Phase coordinates from conventional KDSI may deviate from actual phases as shown in Fig. 1(a), which can be greatly improved in CPF-KDSI [Fig. 1(b)].

3. Result

3.1. System setup

The schematic of the SS-OCT system used for experiment is shown in Fig. 2. A wavelength swept laser (HSL-2000, Santec) with 100nm bandwidth at central wavelength 1.3$\mu$m is adopted and its average output power is 7.9mW. The wavelength swept laser provides a repetition rate up to 100kHz. Spectral data and MZI (INT-MZI-1300, Thorlabs) calibration signal are acquired by a two-channel high-speed digitizer (PCI-5122, National Instrument Inc). A focus-extended fiber probe,²⁷ integrated at the sample arm, scans transversely over samples in cooperation with a 2D linear precision translation stages.

Figure 3 demonstrates the unwrapped phase of the MZI calibration signal and its differentials. It can be observed in Fig. 3(a) that unwrapped phase in start region marked by green box apparently deviates from the constrained polynomial fit curve. In contrast, unwrapped phase in central region, in line with expectation, matches well with the fitting curve. Quantitative analysis is performed by discrete differentials of unwrapped phase [Fig. 3(b)]. Phase in green box labeled region varies at range of 2.1 rad in comparison with only 0.3 rad in the central region marked by orange box.

Once the linear-in-$k$ coordinates are generated, the spline interpolation can be applied to resample the raw spectral data in conventional KDSI or CPF-KDSI. System’s axial PSFs obtained with direct FFT, NDFT, conventional KDSI and CPF-KDSI methods are shown (the fourth-order dispersion compensation is applied during post-processing²⁸) in Fig. 4. All PSFs are normalized at the maximum peak intensity with NDFT. Sliver-coated mirror with 95% reflectance placed at the linear precision translation stage is adopted as the sample for the calibration. Direct FFT causes the broadening of the peak of axial PSF and increases its wings, revealing that the resampling is prerequisite before FFT. Compared with conventional KDSI, CPF-KDSI is capable of achieving narrow PSF and peak intensity close to that with NDFT. One should note that although NDFT can achieve precise results, it is barely capable of real-time imaging for its heavy computational burden. Table 1 specifies the full width at half maximum (FWHM) and normalized intensity of the PSF processed with the four different methods. FWHM of the PSF, or axial resolution at this depth position, is 9.5$\mu$m after CPF-KDSI and FFT, in contrast to 10.6$\mu$m if conventional KDSI is applied. There is also a normalized improvement of 0.1 in peak intensity with CPF-KDSI contrast to that with convention KDSI.

To further investigate the axial resolution and sensitivity enhancement with CPF-KDSI, a plane mirror, as the sample, is located at different depths. Figure 5(a) shows the sensitivity roll-off of the SS-OCT system with conventional KDSI (green curve) and CPF-KDSI (orange curve). Before FFT, the spectrum is zero-padding 64 times to make sure that there are enough sampling points in the depth profiles to precisely determine the axial resolution. From the sensitivity fall-off fitting curve, CPF-KDSI can obtain an improvement of sensitivity over all depths and realize a slight improvement on sensitivity fall-off. Specifically, the 6 dB fall-off of sensitivity with CPF-KDSI is at 2.76mm compared to 2.73mm with conventional KDSI. Quantitative enhancement of the axial resolution and sensitivity with CPF-KDSI is shown in Fig. 5(b). More than 1$\mu$m enhancement of the axial resolution is confirmed over all depths. Furthermore, greater improvements on the axial resolution and sensitivity at deeper position are observed with CPF-KDSI.

3.2. Imaging of microspheres-agar sample and lemon

Homemade polystyrene microspheres-agar sample and a fresh lemon are imaged to interrogate the advantage of the proposed CPF-KDSI. Polystyrene microspheres-agar sample is fabricated by dropping the polystyrene microspheres liquor (nominal diameter 5$\mu$m, concentration 2.5%) to heated gelatinous agar. After solidification, polystyrene microspheres are fixed in agar without noisy motion. Reconstructed OCT images with different methods are shown in Fig. 6. 20 frames are averaged to reduce the speckle noise and enhance SNR.²⁹ Image is severely blurred after direct FFT, as the air-agar surface and boundary of single microsphere are widen drastically in axial direction [Fig. 6(a)]. After resampling with conventional KDSI or CPF-KDSI, the surface and microspheres’ boundary within the focal depth of the probe are much clearer and sharper [Figs. 6(b) and 6(c)]. Diameters of microspheres are smaller than the system’s axial resolution, so sizes of microspheres in the reconstructed images directly reflect the image resolution. By rescaling two areas labeled by green and orange boxes, respectively, within the focal depth region, the microspheres marked by red arrows obviously demonstrate that the CPF-KDSI method results in the minimum size of the microspheres [Fig. 6(d)].

Reconstructed OCT images of a fresh lemon are shown in Fig. 7. The surface is tilted to avoid the strong specular reflection from saturation of the detector and inducing of harmonic noise in the reconstructed images. As is presented in Fig. 7(a), the cell wall of the lemon is obscured under direct FFT. By applying conventional KDSI or CPF-KDSI, the cell wall within focal depth is sharp and of high SNR [Figs. 7(b) and 7(c)]. Furthermore, depth profiles at two positions labeled by green and orange lines reveal that the reconstructed OCT image with CPF-KDSI has better resolution and peak intensity than that with conventional KDSI [Figs. 7(d) and 7(e)].

4. Discussion and Conclusion

In FD-OCT, phase retrieved from the spectrum can be used to obtain an interpolation coordinates for resampling due to its linear relation with wavenumber. However, noise consisted in spectral signal limits the precision of the retrieved phase and therefore degrades the performance of the reconstructed image. A large amount of studies on denoising have been proposed hitherto. Yan surpresses the noise presented in spectral phase with averaging multiple signals.³⁰ Zavareh uses the Kalman approach to filter the noise and generate linear-in-$k$ sampling data in real-time.³¹ All these work demonstrate that denoising is effective to resample spectrum precisely with methods of spectral phase extraction in FD-OCT. In this paper, we use point-constraint to minimize the influence of hazardous noise that is pronounced at the beginning and the end of each sampling.

It is noteworthy that MZI calibration signal is of high SNR, thus phase error induced by noise is not obvious. For this reason, the proposed CPF-KDSI method surpasses conventional KDSI method in a minor improvement in performance of the reconstructed OCT images. This analysis is also testified by the PSFs at different positions. In a deeper position where the PSF has lower SNR, CPF-KDSI can provide a greater enhancement in the axial resolution and sensitivity. In addition, one should note that the constrained polynomial fit must be conducted at each A-scan. Therefore, spectrum resampling with CPF-KDSI consumes more time. Specifically, with our implemented computer (CPU is Intel i5-4570, no GPU accelerating), the time taken for reconstruction of one OCT image ($3800\times 2300$) is 4.8s with CPF-KDSI and 2.6s with conventional KDSI.

SD-OCT is another implementation of FD-OCT, and has no requirement for real-time k-calibration. As same as in SS-OCT, phase extracted from the spectrum can be used for spectrometer calibration in SD-OCT.³² Once the uniform phase coordinates is generated in the process of spectrometer calibration, all spectral data can be resampled with this coordinates. This means that the time costs for both conventional KDSI and CPF-KDSI are almost identical in SD-OCT. The proposed method is thus applicable to SD-OCT.

In conclusion, we propose a $k$-domain spline interpolation method for spectral resampling in SS-OCT. In the presence of phase errors caused by noise and Hilbert transform, unwrapped phases of MZI calibration signal in marginal region derivate from the actual phases. Phase coordinates for interpolation can be obtained much more precisely by polynomial fit with constraint of passing through the central sample point. A comparison of axial PSFs at different depths demonstrate that CPF-KDSI acquires better axial resolution and peak response than conventional KDSI. Reconstructed OCT images of a microspheres-agar sample and a fresh lemon also testify that sharper edges and improved resolution are achieved by the proposed CPF-KDSI method.

Conflicts of Interest

The authors have no relevant conflicts of interest to disclose.

Acknowledgments

The authors acknowledge funding from National Key Research and Development Program of China (2017YFA0700501), National Natural Science Foundation of China (62035011, 11974310, 31927801, 61905214), Natural Science Foundation of Zhejiang Province (LR20F050001), Fundamental Research Funds for the Central Universities.