1. Introduction

Since its first introduction in 1991 [1], optical coherence tomography (OCT) has been developed towards an effective and widely used method for high resolution non-invasive, non-contact and label-free imaging of tissues, fluids, and other scattering media. Especially the introduction of Fourier domain OCT (FD-OCT) with its superior signal-to-noise ratio (SNR) has accelerated the application of OCT imaging as it enables real-time imaging with high resolution and high frame-rate. Currently, it is an established modality in medical imaging [2], and is also applied in other fields such as material science [3,4], artwork examination [5,6] and plant imaging [7].

For imaging small details, the spatial resolution of the imaging system is one of the most crucial aspects. Improving resolution enables imaging of smaller sample structures and more accurate quantification of sample morphology. In OCT, the resolution in the lateral and the axial direction are decoupled. While the lateral resolution depends on the focusing optics and the wavelength, the axial resolution is determined by the temporal coherence of the light source, which is inversely proportional to the source bandwidth. Therefore, improving the lateral and axial resolution require a different approach and are thus often treated separately. In the current work, we focus on improvement of the axial resolution.

In FD-OCT, the sample reflectivity profile along a scan line is usually estimated from the inverse discrete Fourier transform (DFT) of the spectrum, i.e., the interference spectrum and the reflectivity profile form a Fourier transform pair. This limits the (amplitude based) axial resolution to the coherence length of the source, which, for a Gaussian spectrum depends on the center wavelength $\lambda _c$ and source bandwidth $\Delta \lambda$ as $l_c=\frac {2\ln {2}}{\pi }\frac {\lambda _c^{2}}{\Delta \lambda }$. Therefore, OCT setups have been developed using light sources with an ultra-broad bandwidth of up to 200 nm operating at wavelengths typically around 800 nm. Using these approaches axial resolutions up to about a single µm have been achieved [8,9]. However, light sources with an ultra-broad bandwidth are not only difficult to build, their implementation also complicates the optical design of the OCT system. For swept-source OCT (SS-OCT), the source bandwidths are limited due to the limited bandwidth of gain materials and complexity of combining different gain materials [10]. For large bandwidth spectral domain OCT (SD-OCT) the spectrometer performance is the bottleneck, giving a significant intensity roll-off even when properly designed. Moreover, for any ultra-large bandwidth OCT system chromatic aberrations, sample induced dispersion, and wavelength-dependent scattering properties can be hard to compensate for. Another strategy to improve the axial resolution is to use short wavelengths reaching into the visible light range [11]. However, in the blue to green wavelength light range, available superluminescent diodes (SLDs) have a narrow bandwidth and ultra-high resolution cannot be achieved. Therefore, there is a need to overcome this traditional resolution limitation and provide a high axial resolution with a limited source bandwidth.

Early computational methods for improving the axial resolution of OCT were based on classic deconvolution techniques [12,13], with the Wiener filtering and the Lucy-Richardson method being the most notable examples, [14–17]. However, the improvement in the axial resolution achieved by those methods is either limited or obtained at the cost of side effects such as side-lobe or noise enhancement. A significant improvement in axial resolution was obtained with modulated deconvolution [18]. However, this method requires careful calibration of multiple kernels.

High-resolution spectral estimation (SE) techniques, that originate from the field of radar and telecommunication signal processing [19–23] have been used to improve the axial resolution in OCT by using reconstruction methods alternative to the discrete Fourier transform. The frequencies of the interferograms that are measured in FD-OCT correspond to depths of reflection. SE techniques estimate these frequencies with much higher precision than what can be obtained with the conventional DFT, thus effectively improving the axial resolution. Moreover, resolution improvement goes together with excellent side-lobe suppression, yielding high-quality images. Spectral Estimation OCT (SE-OCT) has been able to obtain an improvement in axial resolution with a factor up to 4.7 [21].

Although SE-OCT has shown promising results, there are at least three major limitations that obstruct its usage in the mentioned application areas: the large computational costs, the dependency on user-optimized parameters, and the inaccuracy in reconstructed intensity. Widescale application of SE-OCT techniques would require close to real-time B-scan processing, something that is currently not achieved. Based on reported running time and on the timing of shared code, typical computation times for a single B-scan range from several minutes [21,22] to even hours [23]. Parametric methods only yield accurate results when the model sufficiently fits the data and parameters are well chosen. Optimal parameters often depend on the OCT setup and imaged sample, thus requiring tedious parameters optimization for each new application [19,22]. Even for optimized parameters, misfits between the model and data can cause image artifacts in the results. For example, with the auto-regressive (AR) method, as proposed by Liu et al. [21], the estimated signal intensity is fluctuating and not proportional to the true intensity of back-scattered light, and spurious peaks appear in the images when a high model order is chosen. Similar problems are observed with the maximum entropy method [19] and the Prony method [20]. Other applied SE methods assume a high level of sparsity, thus only yielding accurate results for relative sparse samples [20,23].

Non-parametric methods based on data-adaptive filter-banks, such as Capon [24], Amplitude and Phase EStimation (APES) [25] and the more recent Iterative Adaptive Approach (IAA) [26] form an attractive alternative for parametric spectral estimation methods [27]. These non-parametric methods provide spectral estimations with high resolution and strong suppression of side-lobes. Motivated by the recent advances in non-parametric high-resolution spectral analysis, we propose the use of IAA [26] for the improvement of the axial resolution and image quality in the SD-OCT imaging while avoiding the mentioned limitations of computation costs, dependency on model parameters, and image artifacts. As a non-parametric method, IAA does not require the data to be described by a parametric model, so it does not suffer from a model-data misfit. The computation cost is strongly reduced by using our earlier developed computationally efficient implementation of IAA [28,29]. Additionally, we further reduce computation cost by implementing this algorithm recursively in the lateral direction [30] and limiting the axial reconstruction range to the region of interest, yielding sub-second reconstruction times for a B-scan. With this computational efficient IAA, we aim to make SE-OCT applicable to improve axial resolution in the many application areas where physical accuracy and reconstruction speed are necessary.

2. Theory

2.1 Fourier-domain OCT signal model

The OCT interference signal, i.e., the signal with the constant reference intensity (DC) term, auto-correlation terms, and constant pre-factors omitted, is in Fourier-domain OCT described as [31]

 

Fig. 1. Overview of OCT reconstruction with (a) the conventional DFT-based approach and (b) the iterative adaptive approach (IAA).

Download Full Size | PDF

To set-up the OCT imaging process in the discrete domain we consider the interference signal measured at $N$ discrete wavenumbers $k_n$, where $N$ is normally the number of pixels on the spectrometer camera. Spectrally normalizing this discrete interference signal results in an $N\times 1$ data vector $\mathbf {y}$, with elements $y_n=I(k_n)/S_0(k_n)$ [20,21]. The depth $z$ can be discretized to $M$ uniformly spaced depth locations as

For negligible or uncorrelated (white) noise $\mathbf {\eta }$, $a(m)$ can be estimated as the solution of a least squares data fitting problem of the form [27]

2.2 Iterative adaptive approach method

The iterative adaptive approach (IAA) is a high-resolution, non-parametric spectral estimation method [26,32], proposed in the area of radar signal processing. Figure 1(b) gives a schematic overview of OCT reconstruction with IAA. The reflectivity $a(m)$ is estimated with a weighted least square reformulation of Eq. (5)

As shown in Fig. 1, $\textbf {Q}_m^{-1}$ is initialized as the identity matrix $\textbf {I}_N$, which gives the same result as the DFT with zero-padding. Subsequently, the $N\times N$ matrix $\textbf {Q}_m$ is calculated based on the estimate of $a(m)$ as

The matrix $\textbf {Q}_m$ has a strong weight for signals from depths that have a high amplitude ($\left |a(m)\right |^{2}$), except for the estimated depth $m$, whose contribution is subtracted in Eq. (8). Weighting with the inverse of $\textbf {Q}_m$ thus suppresses the contribution from high-intensity peaks at locations not identical to the estimated signal at location $m$. Thus, side-lobes and edges of the main lobe are suppressed, while the signal that belongs precisely to the location $m$ is passed nearly undisturbed. This results in narrow peaks with limited side-lobes and thus in a high-resolution OCT reflectivity profile $a(m)$. Both $a(m)$ and $\textbf {Q}_m^{-1}$ are further refined during subsequent iterations [26].

Use of $\textbf {Q}_m$ in Eq. (7) results in a considerable computational overload because $\textbf {Q}_m^{-1}$ has to be calculated for every $m$. Fortunately this can be avoided by the introduction of

The noise covariance matrix $\boldsymbol {\Sigma }$ that appears in Eq. (9) has to be estimated from the given data. A common assumption is that $\boldsymbol {\Sigma }$ is a diagonal matrix with (non-negative) elements $\sigma ^{2}(n), n=1,\ldots, N$. These can easily be estimated, in a scheme similar to Eq. (10), as

The IAA estimator is obtained by iterating Eq. (9) and Eq. (10) until convergence. About 10 iterations are usually required for convergence, see also [26,32,33]. The algorithm is summarized in Appendix A.

2.3 Brute force IAA

When implementing IAA in a brute force way, i.e., by direct use of Eqs. (9), 10, 11 and 12, without taking into account the underlying structure of the pertinent matrices and variables, the overall computational complexity denoted with $\mathcal {C}$ is approximately given by

2.4 Fast IAA

The computational complexity can be significantly reduced by taking into account the particular properties of the variables resulting from the Fourier vectors (Eq. (4)) [28,29,34]. As we assume the noise variances to be equal, i.e., $\boldsymbol {\Sigma }=\sigma ^{2}\textbf {I}$, the interference covariance matrix $\textbf {R}$ is a Hermitian Toeplitz matrix, which instead of using Eq. (9), can be efficiently estimated by means of Toeplitz to Circulant matrix embedding in ${ \mathcal {O}}\left (M\textrm {log}_2(M)\right )$ operations [28]. Taking into account the special structure (the low displacement rank) of a Toeplitz matrix, $\textbf {R}^{-1}$ can be implicitly estimated by means of the Gohberg-Semencul (GS) factorization and the use of the Levinson-Durbin algorithm [27], briefly explained in Supplement 1. Using the GS factorization of $\textbf {R}^{-1}$, the elements in Eq. (10) and the noise variance elements $\sigma ^{2}(n)$, Eq. (11), are computed using the FFT as the computational engine, as it is shortly explained in Supplement 1. This gives an overall computational complexity of approximately [28,29,34]

2.5 Recursive scheme for B-scan processing (RFIAA)

B-scan OCT imaging is performed by processing $N_B$ consecutive A-scans as columns in an image matrix. The columns are usually processed independently to produce the corresponding sequence of depth profiles, eventually combined to an OCT B-scan image. Obviously, the computational cost for processing a B-scan using the IAA method for independent A-line processing is given by

In the most efficient algorithm, RIAA is implemented based on the FIAA method, resulting in the recursive fast IAA (RFIAA) method. The recursive scheme is made compatible for parallel processing by dividing the B-scan into data chunks according to the number of available CPU cores. For each chunk, the first column is initialized with $q_i$ iterations, and subsequent columns are iterated $q_{rci}$ times after the proposed initialization close to convergence. As the number of cores is usually much smaller than $N_B$, the computation complexity is not significantly affected and still follows Eq. (15). However, this implementation makes RFIAA applicable with parallel processing with multiple cores which results in significant computation time reduction.

2.6 Reduction of reconstruction range

The standard FD-OCT reconstruction implies a 50% redundancy in the estimated depth profiles as the mirror image (with negative $z$-values) is usually rejected, keeping only the values for the positive depth range. Moreover, due to the presence of auto-correlation artifacts and limited imaging depth, the depth range that contains useful sample information is usually even smaller.

We increase the efficiency of the OCT reconstruction with IAA by restricting the reconstruction to a useful depth range $\Delta z$. This reduces the reconstruction grid length $M$ with a factor $1/R=\Delta z/(2z_{\textrm {max}})\leq 0.5$. Limiting the reconstruction depth range also allows for reduction of the input interference spectrum length $N$ such that it only includes the frequencies that contribute to the useful range. This is implemented by shifting the useful depth range to the region centered at $z=0$ using the modulation property of the DFT, followed by low-pass filtering and down-sampling at the rate $1/R_s$, where $R_s=\left \lfloor R\right \rfloor$ is the integer value operator. This reduction is most efficient when $R$ is an integer that equals a power of 2. The computational gain due to range reduction is a factor $\mathcal {O}\left (R^{2}\right )$ or $\mathcal {O}\left (R\left (1+\log _2\left (R/M\right )\right )\right )$, depending on whether the first or the last term in Eq. (14) is predominant.

3. Methods

3.1 Experimental setup

The measurements were performed with a commercial spectral-domain OCT system (Ganymede-II-HR, Thorlabs, Germany). The light source consists of two coupled super-luminescent diodes that are combined in a spectrum with a central wavelength of 900 nm and a full width half maximum (FWHM) of 120 nm. The FWHM of the axial intensity point spread function (PSF) was experimentally measured to be 3.0 µm in air, which corresponds to 2.2 µm in tissue. The spectrometer has 2048 pixels covering a 220 nm bandwidth, giving an axial imaging range of 1.87 mm. The intensity lateral PSF was experimentally determined with a knife-edge step response to have a FWHM of 6 µm in focus. The spectral raw data was acquired using ThorImage software (version 5.2.0) and further processed in MATLAB (R2020a).

3.2 OCT sample imaging

Experimental OCT data was obtained from a wedge phantom, a layered interface phantom, an onion sample, and finger skin tissue.

The wedge phantom represents a sparse object that is ideal to study the ability to resolve two closely separated reflectors. The air wedge was constructed by placing a coverslip on a microscope glass; it was tilted by placing a piece of tape between them at one end. A neutral density filter (NE20A-B, Thorlabs, Germany) was placed between the lens and the wedge phantom to circumvent detector saturation, and the wedge was placed in a depth region without auto-correlation artifacts. The position of the two wedge interfaces was determined by a two-peak Gaussian fit on the FBW-DFT A-scans in the region where the interfaces were well resolved. The interface location was extrapolated using a quadratic fit to the region where the two interfaces were not well resolved. This extrapolation served as the ground truth to determine the separation of the interfaces.

The layered interface phantom was used to study intensity fluctuations at different SNRs. It consisted of two coverslips on a microscope glass. The coverslips were spaced by a layer of tape at their edges and the space was filled with ethanol to create both air-glass and ethanol-glass interfaces with different intensities and SNRs. The phantom was placed out of focus to avoid detector saturation.

To test IAA on medium sparse objects, a slice of an onion was imaged. Onion cells are typically 50-100 µm high and have clear boundaries that consist of two closely separated interfaces, which are typically 6-10 µm separated [18]. Thus, these layers can only be resolved with a high axial resolution. First, a cross-sectional image was taken that contains 1024 scan lines over a range of 3 mm, giving a lateral spacing of 2.9 µm between the scan lines. Second, a volume data set was obtained of $512\times 512$ scan lines, covering $2.5\times 2.5$ mm$^{2}$. With the volume scan, 8 spectra were acquired per lateral position, which were averaged before further processing to increase SNR.

IAA was tested on a non-sparse sample by imaging skin tissue from the fingertip. The dense skin tissue has a low level of sparsity and is a good sample to evaluate the performance of the methods on reconstruction quality and contrast in the context of medical imaging. The B-scan comprises 2048 lines over a lateral range of 8.1 mm, giving a lateral spacing of 3.9 µm between subsequent scan lines. The contrast to noise ratio (CNR) was calculated as

Finally, contrast and speckle statistics was quantified using an Intralipid suspension (Fresenius-Kabi) diluted to 2.5 weight%. A droplet of the suspension was placed under a microscope coverslip to create a flat top surface, which was aligned to be horizontal in lateral scan direction. A B-scan of 1024 scan-lines was recorded over a lateral range of 3 mm. The depth dependent intensity due to confocal PSF, roll-off and attenuation [35] was compensated for by dividing the OCT intensity of the Intralipid region by the laterally averaged intensity of this region, thus creating a homogeneous scattering region.

3.3 OCT data simulations

OCT data was simulated to study the resolution, intensity preservation, and contrast in RFIAA reconstruction. The simulated spectra were based on a 1D OCT model [31]

 

Fig. 2. (a) The original source spectrum and the reshaped spectra for each of the four methods. (b) A typical zoom in A-scan obtained with RFIAA for different $q_{rci}$.

Download Full Size | PDF

For the resolution study, a wedge was simulated as two reflectors with an equal amplitude and a spacing $z_2-z_1$ ranging from 0 to 29 µm over 1024 A-scans. The simulated object for studying intensity preservation consisted of 8 horizontal interfaces with halved amplitude (-6 dB intensity drop) for each subsequent interface. For studying contrast and speckle statistics, 3 speckle regions with different mean intensity were simulated. The spectrum from a speckle region was simulated by taking 2048 reflectors with equal amplitude $a_j$ located at random depths $z_j$ (uniform probability distribution) over a 150 µm depth range. The signal from these sub-resolution reflectors added up coherently, creating a 1D speckle pattern. The amplitudes $a_j$ were halved for each subsequent region. The position of the reflectors in subsequent A-scans was independent, so the effect of finite lateral resolution is not taken into account.

3.4 OCT data processing

Pre-processing of the spectral data consisted of three steps: subtracting the reference arm OCT spectrum, interpolating the data on a grid that is linear in the wave-number domain, and multiplying with a phase vector to compensate for dispersion mismatch. The dispersion was corrected for with a fourth-order polynomial, whose coefficients are determined from a reference measurement with a single mirror as reflector [36]. After the generic pre-processing, the data were further processed according to four different methods:

3.4.1 Full bandwidth-DFT (FBW-DFT)

This is the DFT reconstruction of the full bandwidth signal and acts as a ground-truth image. The interference spectra were Gaussian reshaped to suppress side-lobes [37], as shown in Fig. 2 (a). To avoid noise amplification, the edges the spectrum were not reshaped but followed the original source spectrum. The FWHM of the corresponding bandwidth-limited axial PSF is 2.7 µm. The OCT depth information was obtained by taking the inverse DFT of the reshaped spectral data, with zero padding to obtain a grid spacing that is equal to that of the RFIAA reconstruction.

3.4.2 Partial bandwidth-DFT (PBW-DFT)

This is the DFT reconstruction of a partial bandwidth interference spectrum that is obtained by reshaping the spectrum with a narrow Gaussian window. The spectral interval was chosen around the high-intensity peak on the left side of the original spectrum, as indicated in Fig. 2. This ensured that the majority of the signal comes from a single SLED, thus avoiding non-uniform modulation in the region where the spectra of different SLEDs overlap. Non-uniform modulation can be caused by polarization mismatch between the SLEDs, thereby degrading the axial resolution [38] and potentially affecting the spectral estimation result. As low-intensity spectrum tails have insufficient SNR for spectral estimation [21] but help in side-lobe reduction for DFT reconstruction, we chose the Gaussian window for PBW-DFT to have a relative intensity of 10% at the edges of the window that is used for the spectral estimation methods and extend it to outside this domain, as shown in Fig. 2. The FWHM of the corresponding bandwidth-limited axial PSF is 8.3 µm. The reshaped spectral data was zero-padded before taking the inverse DFT to obtain a grid spacing equal to that of the RFIAA reconstruction.

3.4.3 Auto-regression filter (AR171)

This is the auto-regressive reconstruction on a spectrally normalized truncated part of the bandwidth that covers a quarter of the spectrometer bandwidth, identical to the bandwidth of the PBW-DFT except for the low intensity edges, as indicated in Fig. 2. The AR parameters were estimated using the modified covariance method and the model order was set to 171, corresponding to a third of the data pixel length as proposed by Liu et al. [21]. The DFT-based axial PSF of the corresponding rectangular source spectrum has a FWHM of 6.5 µm.

3.4.4 Recursive fast IAA (RFIAA)

This is the proposed RFIAA spectral estimation OCT on the truncated part of the total bandwidth (same part as for AR171). The spectral data is spectrally normalized before using it in the RFIAA algorithm, resulting in a DFT-based axial PSF with a 6.5 µm FWHM.

The number of iterations $q_i$ of the first scan line was set to 10, which was sufficient to obtain convergence. The value of $q_{rci}$ was chosen based on a test with a reconstructed A-line for different $q_{rci}$ of the onion sample, shown in Fig. 2(b). A value of $q_{rci}=2$ already gives significant side-lobe reduction (1-3) and deeper valleys (4) over $q_{rci}=1$. The further effect of increasing $q_{rci}$ beyond $q_{rci}=2$ is limited. Hence, we chose $q_{rci}=2$. The RFIAA upsampling factor $M/N$ (for full range reconstruction) was chosen to be 64 for the wedge and simulation images to achieve an optimal evaluation of the performance, especially for high SNRs that give a high resolution. For the onion and skin sample, $M/N=16$ was chosen for fast computation and good image quality. The reduced reconstruction depth range for all objects covered half the positive $z$-range (downsampling factor $R=4$) centered around region of interest.

For all four methods, the intensity of the OCT signal is used for further analysis. For visualization, the intensity is log compressed and the lower limit of the dynamic range is determined by fitting a Rayleigh distribution on the histogram of amplitudes of a custom selected noise region in the image, as described by Steiner et al. [39] (using $\tau =0.95$).

All data processing was implemented as MATLAB scripts (version R2020a) using its intrinsic functions. FIAA implementation is taking advantage of the levinson.m intrinsic function, which is actually implemented as a mex file bundled in the MATLAB programming environment. No further attempt for speeding up the running time using MATLAB’s compiler capabilities nor any direct C-code programming has been applied. The AR modified covariance estimator, which is required in the approach proposed in [21], is implemented using armcov.m intrinsic function. This is a plain, brute force implementation of the AR modified covariance estimator, offering some intrinsic parallel processing capabilities due to the use of the available BLAS3 routines. All code was executed on a Dell Precision 5820 with an Intel Xeon W-2223 CPU and 32 GB RAM. The B-scan was cut in 4 equal-sized sub-images for parallel processing on the 4 CPU cores, as discussed in 2.5. AR reconstruction was applied in parallel on scan line level.

4. Results

4.1 Sparse wedge object

Figure 3 shows the imaging results of the wedge phantom obtained with the different methods. The two interfaces of the wedge are clearly visible (Fig. 3(a-d)). The FBW-DFT image shows a beating pattern at the location where the interfaces approach each other. This effect becomes a lot more pronounced when the bandwidth is reduced with PBW-DFT. This beating pattern is caused by coherent addition of the reflection of both interfaces, and further influenced by spectral leakage in DFT reconstruction [40]. For PBW-DFT, the beating pattern strongly reduces the resolvability of the two interfaces and is an indication of the lower axial resolution. Using the same bandwidth as with PBW-DFT, the two spectral estimation methods reconstruct the two wedge interfaces as narrow lines and largely eliminate the beating pattern. AR171 gives the most narrow representation of the interface, which appears relatively dark due to strong variation of peak intensities caused by the AR reconstruction. This variation required the dynamic range of the image color scale to be extended to 55 dB. Where the interfaces approach each other, vertical stripe artifacts appear and the interfaces become ill-defined. Though RFIAA gives a slightly less narrow representation of the interface than AR171, the intensity of the interfaces are more constant. Still the improvement with respect to PBW-DFT is very clear and the interfaces are imaged narrower than FBW-DFT that uses a four times larger bandwidth. Towards the left of the image, the interfaces melt together in a beating pattern, and some stripe artifacts are visible.

The cross sections in Fig. 3(e) further illustrate the described effects. PBW-DFT gives a broad peak from the merged interfaces. AR171 gives narrow peaks, but the intensity of the first interface is less than 10% of the second interface which does not resemble the true relative intensity as obtained with FBW-DFT. This effect is clearly a disadvantage of AR reconstruction [21] and will be explored further in the next section. RFIAA gives peaks that are about the same height and a bit narrower than FBW-DFT.

 

Fig. 3. (a-d) Experimental OCT images of the wedge phantom. The SNR is 32 dB, and the dynamic range is 15 dB (a,b,d) or 55 dB (c). Scale bars indicate 10 µm vertical and 100 µm horizontal. The arrows indicate the widest separation where the surfaces are not resolved. (e) Cross section of the OCT intensity at the white dotted line in (a-d). (f) The resolution as a function of SNR based on simulations (solid line). The experimental values are indicated with the star.

Download Full Size | PDF

The OCT axial resolution was quantified as the widest spacing at the edge where the interfaces first merge for two successive A-scans. The two interfaces are considered to merge when the valley between the peaks is less than a factor 0.5 (3 dB) below the lowest peak intensity. The arrows in Fig. 3(a-d) indicate the resolution limit for the shown images, at a spacing of 4.9 µm, 14.4 µm, 4.2 µm and 5.5 µm for FBW-DFT, PBW-DFT, AR171, and RFIAA respectively. RFIAA thus gives a 2.6 times better resolution than PBW-DFT that uses the same spectral bandwidth. Figure 3(f) shows the resolution as a function of SNR, based on simulated OCT data of a wedge. The resolution for DFT reconstruction is independent of SNR and thus follows a horizontal line. For spectral estimation methods, the resolution improves with increasing SNR [21], which is visible in the decreasing curves for AR171 and RFIAA. For an SNR above 30 dB, AR171 outperforms RFIAA and achieves resolutions up to 0.2 µm for SNR > 75 dB where RFIAA obtains a resolution of 1 µm. However, in practice, few OCT measurements are performed at such a high SNR. For SNR values below 30 dB, RFIAA outperforms AR171, showing that it is better able to handle data with a medium to low SNR. Over the full studied SNR range, RFIAA gave a better resolution than PBW-DFT. The resolutions from the experimental data are indicated in the graph with a star and correspond well to the simulated results.

4.2 OCT intensity reconstruction

A major drawback of AR171 is that reconstructed intensity does not correspond with the true intensity of the reflected light. In this section, we further investigate this behavior and compare it to the performance of RFIAA. Figure 4(a) shows the image of an OCT simulation of 8 interfaces whose intensities drop with 6 dB with every subsequent interface. The decreasing intensity is clearly visible in the DFT reconstructed regions of the image. Figure 4(b) shows part of an A-line at the dashed lines in (a) for the different methods. For FBW-DFT, PBW-DFT, and RFIAA, the peaks 1 to 6 have monotone decreasing intensities corresponding to the simulated reflection intensities. In contrast, the peaks of AR171 do not follow this pattern and give fluctuating intensities that do not correspond to the real reflection intensity; e.g., peak 6 has a 10 dB higher intensity than peak 1, while it should be 30 dB lower than peak 1.

This effect is also clear when we take all the 1024 lateral positions into account, as shown in Fig. 4(c). Here, the mean intensity of each interface is indicated with a diamond, together with error bars that mark the edges of the middle 95% intensity values, i.e., the $2.5^{\textrm {th}}$ and $97.5^{\textrm {th}}$ percentile. These edges correspond to $2\sigma$ below and above the mean for normally distributed values. On the horizontal axis is the average SNR of the corresponding interface in the PBW-DFT image. The mean intensities for FBW-DFT, PBW-DFT, and RFIAA lie on a straight line with a slope of 1 and have small error bars. AR171 gives much lower mean intensities, due to a few high-intensity peaks on which the image is normalized. Moreover, the mean intensities do not follow a straight line, but flatten out for high SNR above 30 dB, showing that even the average intensity over 1024 lines does not follow the true reflection intensity. Most importantly, the intensity varies with 20 dB to 30 dB, clearly showing that AR reconstruction does not allow for quantifying relative reflection intensities of features within a sample.

 

Fig. 4. Results on OCT reflection intensity preservation. (a) OCT reconstruction of a simulated object with interfaces whose reflection intensity drops with 6 dB with every subsequent interface. The different columns indicate the results of the four methods. The scale bar is 100 µm. (b) Part of an A-scan at the dashed lines in (a), normalized to the maximum value of the A-scan. (c) The mean OCT peak intensity of the interfaces. Error bars mark the range of the middle 95% intensity values. The SNR is determined from the PBW-DFT image. (d) Full range of the middle 95% intensity values as a function of SNR of the PBW-DFT image.

Download Full Size | PDF

The width of the middle 95% intensities range (95% width) as function of SNR is plotted in Fig. 4 (d). PBW-DFT corresponds well with the theoretical value, which is based on the standard deviation of peak intensity and noise [41], and is indicated with the black dashed line. This shows that all the fluctuations are caused by the noise and not by reconstruction inaccuracies. FBW-DFT is below this theoretical value because with a larger bandwidth it has more power and thus a higher SNR than PBW-DFT, whose SNR is on the horizontal axis. The 95% width with RFIAA is less than 0.5 dB above the theoretical value, a difference that is negligible compared to typical dynamic ranges in OCT images. It is worth noticing that for high SNR, where the theoretical values and DFT-based methods go to 0, RFIAA remains at a plateau of about 0.5 dB. This may indicate that where DFT-based methods reduce the variation in peak intensity, RFIAA uses this extra SNR to increase the resolution while keeping the variation in peak intensity at an acceptable level. For the DFT and RFIAA methods, the 95% width remains below 3 dB for SNRs above 20 dB. This is in sharp contrast with AR171, that has 95% widths between 20 dB and 30 dB, which are outside the vertical range of Fig. 4(d).

The simulation results on intensity reconstruction were complemented by measurements on the layered interface phantom. The four interfaces had SNRs ranging from 18 dB to 32 dB. The mean intensities and 95% widths are indicated with stars in Fig. 4(c-d) and correspond well with the simulations.

In conclusion, RFIAA, in contrast to to AR171, reconstructs true reflection intensity and thus allows for quantitative analysis based on intensity (e.g. measuring optical attenuation).

4.3 Medium sparse sample

As a biological object, we used a slice of onion tissue, which has large cells without internal scattering structure, except for the nucleus, thus having a medium sparsity level. The cell walls are made up of two layers that are typically 5-11 µm (6.5-15 µm optical path length (OPL)) apart, and could thus only be resolved with a high axial resolution.

Figure 5(a-d) shows the OCT images of the slice of onion, clearly revealing the cellular structure. In the FBW-DFT, the two cell wall layers are clearly distinguishable (indicated by the white arrows), though side-lobes around high-intensity reflections (indicated by the green arrow) slightly reduce the image quality. The limited resolution of the PBW-DFT reconstruction causes the double layers to merge into a thick cell wall. Even at the top surface where the separation between the interfaces is larger, they are barely resolved. AR171 and RFIAA obtain a much higher axial resolution with the same bandwidth as PBW-DFT. AR171 gives the most narrow lines, but the earlier described intensity fluctuations sometimes obscure one of the layers (e.g. at the location indicated by the right white arrow). At this position, RFIAA gives a clear image that is similar in resolution to FBW-DFT. RFIAA is designed to improve resolution while suppressing side-lobes, which is visible in the low side lobes around the high-intensity reflections (compared to FBW-DFT).

 

Fig. 5. OCT imaging results of the onion and skin sample. (a-d) are images of the onion, (e-h) are skin images and (i) and (j) are A-scans from respectively the onion and skin image at the dotted lines. The images are visualized in dB scale with the dynamic range lower limit determined from noise statistics. Scale bars indicate 150 µm. 3D data of the onion and skin sample are shown in Visualization 1 and Visualization 2.

Download Full Size | PDF

The A-scan in Fig. 5(i) gives further insight into the performance of the studied methods. At interface 1, only PBW-DFT does not clearly resolve the two layers. AR171 gives the most narrow peaks, but the second peak has a much lower intensity and is thus barely visible in the image. Although RFIAA gives wider peaks, they are at the correct FBW-DFT intensity and the two layers are clearly resolved. In between the peaks, AR171 gives a smooth valley which is typical for AR spectral estimation. However, as these valleys are not smooth in the lateral direction, vertical stripes are visible as side-lobes in the image, which is enhanced by the necessarily large dynamic range. At interface 2, RFIAA performs well with a deeper valley and less side-lobes than FBW-DFT. AR171 again gives a much lower peak intensity and with a shallow valley, the image does not show clearly distinguished peaks. The distances (in OPL) between the peaks for FBW-DFT, AR171, and RFIAA respectively are estimated at 14.4 µm, 17.1 µm and 14.0 µm for interface 1, and 6.4 µm, 6.0 µm and 7.3 µm for interface 2. This corresponds well with each other and only AR171 on the first interface is more than 1 µm different from FBW-DFT. Concluding, RFIAA improves the resolution significantly over PBW-DFT and gives an image close to FBW-DFT with slightly less side lobes. In contrast to AR171, RFIAA reconstructs reflection intensity in better agreement with FBW-DFT which leads to enhanced image quality.

4.4 Non-sparse skin sample

The second biological object was a skin sample from a fingertip that has a low level of sparsity. Figures 5(e-h) show the OCT images of the skin sample, with the epidermal (ep) and dermal (de) layer, and the helical-shaped sweat duct (inset) clearly visible. AR171 and RFIAA have a higher axial resolution than PBW-DFT, visible in the top surface and axial thickness of the sweat duct. The top left corner of the inset shows a small crack in the top surface that is better resolved with AR171 and RFIAA. Using the red box as signal area and the green box as noise area, the CNRs are 0.76, 0.80, 0.47 and 0.64, for FBW-DFT, PBW-DFT, AR171, and RFIAA, respectively. RFIAA reconstruction gives a higher CNR than AR171, which is also visible from the appearance of the images. Moreover, where RFIAA gives well-developed speckle, AR171 gives some spurious peaks that form narrow curly structures that may be mistaken for sample features.

The A-scan at the white dotted line, Fig. 5(j), shows the peaks corresponding to the skin surface (number 3) and the sweat duct interfaces (indicated by numbers 4,6,7). PBW-DFT gives the broadest peaks, which RFIAA manages to narrow down to a level similar to FBW-DFT (numbers 3 and 4) or even smaller (numbers 6 and 7). AR171 gives the narrowest peaks. However, peak 7 has a significantly lower intensity than the true reflected light intensity, causing it to be less narrow and less distinct from the side lobes. RFIAA has a few high-level side lobes or speckles (number 5) that are still 7 dB below the lowest peak (6).

4.5 SE-OCT CNR and noise statistics

Figure 6 shows experimental and simulation results of homogeneous speckle regions. Figure 6(a) shows the images of the Intralipid suspension, which gives a homogeneous speckle pattern. As speckle size is closely linked with the spatial resolution, RFIAA gives a finer speckle pattern than PBW-DFT. The RFIAA speckle size is slightly coarser than for FBW-DFT, probably because the SNR of a single speckle that consists of a combination of multiple, unaligned, sub-resolution reflectors is relatively low. With AR171, the speckle region appears dark due to some high-intensity peaks that bring the dynamic range down. These peaks in the speckle region also cause a high variance, resulting in a CNR that is 70% below that of PBW-DFT. RFIAA gives a much better contrast with a CNR that is only 15% below that of PBW-DFT. Figure 6(b) shows the speckle amplitude distribution, which follows the predicted Rayleigh distribution for the DFT reconstructed images [42]. The histogram for RFIAA is very close to a Rayleigh distribution, while for AR171 it deviates significantly from this theoretical distribution. Each method gives similar histograms for the amplitude distribution in the noise region.

 

Fig. 6. Experimental and simulation results on contrast and speckle statistics. (a) images of an Intralipid suspension with their CNR. (b) the histogram of the OCT amplitudes in the speckle region with a fit of the Rayleigh distribution (red line) (c) simulation images of 3 regions with different speckle intensities and their CNR. (d) histograms of the OCT amplitude in the noise region with a fit of the Rayleigh distribution (red line). The green dashed lines indicate the lower dynamic range limit for the images of (c). Scale bars indicate 200 µm, where the lateral size of the images in (c) is undefined. The white dashed rectangles indicate the noise and speckle regions for CNR calculation.

Download Full Size | PDF

Figure 6(c) shows the simulated images of three speckle regions with decreasing intensities. FBW-DFT, PBW-DFT, and RFIAA give similar image contrast between the regions, with the size of the speckles being the main difference. The speckle size is similar to that in Fig. 6(a). The CNR of AR171 is 93% to 99% below the CNR of the DFT methods, while the CNR with RFIAA is only 10% below that of the DFT methods. The simulation results on speckle regions of different intensities correspond to the observations from the experimental data. Figure 6(d) shows the histograms of the OCT amplitude in the noise region (top rectangle in Fig. 6(c, FBW-DFT)), which have similar shapes as those for the speckle region in Fig. 6(b). These results indicate that RFIAA outperforms AR171 in preserving noise and speckle statistics, which is useful for, for example, tissue characterization [43] and automatic thresholding [39]. Moreover, they show that as a non-parametric method, RFIAA describes the underlying physics better than AR spectral estimation that assumes an AR model with a fixed order.

4.6 SE-OCT computation time benchmarking

Table 1 shows the computation time for the onion slice sample, averaged over 100 realizations, along with a theoretical count of the complex-valued operations required for processing of a single A-scan by each tested method.

FBW-DFT and PBW-DFT have the same computation time as both are zero-padded to the same length. As expected, DFT reconstruction is the fastest, with a computation time of one to two orders of magnitudes below the other methods. RFIAA with parallel processing on half the useful range ($R=4$, $N=128$, $M=2048$) gives a fast 0.37 s computation time for a 1024-line B-scan. Although this is not yet sufficient for video-rate processing at 20 frames per second, it approaches real-time reconstruction while for post-acquisition viewing the waiting time is negligible. It also allows for fast 3D (see Visualization 1 and Visualization 2) and time-lapse OCT reconstruction within short time scales without the need for supercomputers.

Table 1. Average computation time of a 1024-line B-scan for the studied methods over 100 realizations. The reduced reconstruction range and parallel processing over 4 CPU cores are optional as indicated. AR43 acts on the same reduced range data as RFIAA with reduced range. For the calculation times in this table, $N=512$, $M=8192$ and $R=4$, as used for the onion and skin sample.

View Table

AR171, with 5.23 s computation time, is significantly slower than RFIAA. Reducing the reconstruction range similar to that used for RFIAA will bring the computation time down to 0.43 s for AR43 (AR with an order of 43), but the resulting lower order may hamper the ability to reconstruct complex sample structures [21]. We note that the AR modified covariance estimator could have been implemented using a fast ${\mathcal {O}}\left ( p^{2} \right )$ algorithm [44] instead of the current implementation. However this is neither available in the MATLAB programming environment, nor is it supported as sharing code in [21].

Parallel processing over 4 CPU cores improves the computation speed only by about a factor two for AR171 and both RFIAA implementations. Though the improvement with parallel processing is limited, probably due to overhead with handling large matrices, it still helps to bring the computation time down and allows for fast spectral estimation OCT. These results show that our efficient RFIAA implementation of SE-OCT gives fast OCT reconstruction; it is faster than standard AR spectral estimation implementation, especially for full depth range reconstruction.

5. Discussion

We developed RFIAA, an optimized implementation of IAA, to successfully address three problems that arise in SE-OCT:

In this discussion, we first compare our method with other SE-OCT methods and explore limitations of our method and potential improvements.

Comparing our method quantitatively with other SE-OCT methods is not straightforward. The obtained image quality depends on variables like bandwidth, SNR, imaged sample, and dynamic range of visualization, which differs between publications. In addition, the computation time depends on computer hardware and software implementation. As a way forward, the algorithms should be tested on the same data set, as was done with AR in the current work. To allow for comparison of future work, our data is available in a repository [45].

Still, a cautious and more qualitative comparison with earlier published methods can be made based on the results they report. The iterative re-weighted approach, as proposed by Mousavi et al. [22] is similar to RFIAA, with the main difference that it uses a different weighting matrix, which depends on user-optimized parameters. Moreover, it needs up to 50 iterations, does not adopt a recursive scheme, and uses brute force matrix inversion, leading to A-scan processing times between 1.7 and 18.2 s (29 min to 5 hours for a 1024-line B-scan). Ling et al. [23] used a $\ell 1$-norm minimization, to promote sparsity of the solution, together with a non-weighted minimum least square solver. This method needs a user-chosen Lagrange multiplier and has a large computation load with processing times of 4 s per A-scan (1 h for a 1024-line B-scan) on a high-performance computer. The resolution improvement is similar to what was obtained with RFIAA, though comparison is again not straightforward because they do not mention SNR levels for the wedge experiment and use a different resolution criterion. Both methods allow for faster implementation, but that would need significant extra work and alteration of the algorithm. In conclusion, RFIAA obtains similar results to two recently published state-of-the-art methods without the need for parameter optimization and with orders of magnitude less computation time.

Further improvement in computation time is possible using improved computation hardware with a higher level of parallelization. Data-parallel testing of the proposed RFIAA approach on a 10-core workstation (Intel core i9-9900X) provided significant speedup, giving a fast 0.0607 s computation time for a 1024-line B-scan, in the case where RFIAA with parallel processing on half the useful range ($R=4$, $N=128$, $M=2048$) is considered. These timing results indicate that the proposed RFIAA approach is suitable for parallel implementation on multi-core CPUs, even when fast prototyping, using high-level programming is considered. In terms of a possible future direction, the above results are indicators that a lower-level programming implementation (e.g. C/C++), along with parallel processing pragmas (e.g. OpenMP, OpenACC) deployed on multi-core CPUs and/or accelerators (GPUs) may be the way towards a real-time operating end application.

In this context, it should be highlighted that the most time-consuming part of the proposed RFIAA approach is the linear system solver, which, due to the Toeplitz structure of the underlying linear system, is tackled here using the ${\mathcal {O}}\left ( p^{2} \right )$ Levinson’s algorithm available in MATLAB. However, several other alternative algorithms can be used instead, offering either further computational reduction or being suitable for parallel implementation. Among them, are the so-called superfast ${\mathcal {O}}\left ( p\log _2(p)^{2} \right )$ methods [46–49], the ${\mathcal {O}}\left ( p\log _2(p) \right )$ techniques based on iterative linear solvers such as the conjugate gradient method [50,51], and finally the fully parallelizable Schur-type algorithms, amenable for implementation on massively parallel hardware, [44,47]

SE-OCT works especially well for samples that have a high level of sparsity and clear interfaces with a high SNR. RFIAA is however more robust to high noise levels than AR spectral estimation, as it gives a better resolution for SNR < 25 dB. It also obtains a reliable reconstruction with good contrast of non-sparse regions in a sample, though image quality improvement with respect to DFT based on the same bandwidth is disputable in these areas. However, for a partially sparse sample, RFIAA is able to sharpen the high-SNR features while retaining a good contrast reconstruction of non-sparse areas. As the resolution is SNR dependent, the SNR could be enhanced by spectral averaging, enabling higher resolution at the cost of longer acquisition times.

In this work, we assumed equal variance for all the spectral noise components. This is implemented in Eq. (12), where $\sigma ^{2}(n)$ is averaged and multiplied with the identity matrix to obtain a diagonal with constant values. However, IAA also can use a wavenumber-dependent noise variance $\sigma ^{2}(n)$. The implementation of this would allow for accurately incorporating the higher noise at lower intensity edges of the spectrum while reducing the influence of the higher noise level for these parts. However, implementation of this makes the data-covariance matrix $\textbf {R}$ non-Toeplitz. Consequently, the fast implementation of IAA needs to be adapted leading to an increase in computational load. Thus wavenumber dependent noise variance is most valuable when computation speed is less important than precision.

Contrary to AR spectral estimation, RFIAA obtains both the amplitude and phase of the signal. This would allow for the RFIAA also in the domains of Doppler OCT imaging and sub-resolution phase-resolved OCT motion imaging [52].

6. Conclusions

This paper presented RFIAA, a fast implementation of the non-parametric iterative adaptive approach, for the reconstruction of OCT images with significantly better axial resolution than conventional DFT reconstruction. This SE-OCT method successfully addresses three problems of previously developed SE-OCT, namely the dependency of the reconstruction result on subjective user-set parameters, the occurrence of reconstruction artifacts, and the large computational load. Contrary to AR spectral estimation, the non-parametric RFIAA is consistent in reconstructed intensity, yields a high contrast, and shows less spurious peaks. With a reconstruction time of 0.37 s for a 1024-line B-scan, RFIAA is significantly faster than other SE-OCT presented in literature. This brings SE-OCT a significant step closer to application.