Open Access
Add to BibSonomy
Abstract
Phase-preserving spectral estimation optical coherence tomography (SE-OCT) enables combining axial resolution improvement with computational depth of field (DOF) extension. We show that the combination of SE-OCT with interferometric synthetic aperture microscopy (ISAM) and computational adaptive optics (CAO) results in high 3D resolution over a large depth range for an OCT system with a narrow bandwidth visible light super-luminescent diode (SLD). SE-OCT results in up to five times axial resolution improvement from 8 µm to 1.5 µm. The combination with ISAM gives a sub-micron lateral resolution over a 400 µm axial range, which is at least 16 times the conventional depth of field. CAO can be successfully applied after SE and ISAM and removes residual aberrations, resulting in high quality images. The results show that phase-preserving SE-OCT is sufficiently accurate for coherent post-processing, enabling the use of cost-effective SLDs in the visible light range for high spatial resolution OCT.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
Corrections
6 January 2024: A correction was made to reference 40.
1. Introduction
The spatial resolution in optical coherence tomography (OCT) determines the size of the smallest details in the sample that can be visualized. In OCT, the axial and lateral resolution can be considered to be decoupled [1]. The axial resolution is determined by the coherence length of the light source, which scales with the center wavelength $\lambda _c$ and source bandwidth $\Delta \lambda$ as $\lambda _c^2/\Delta \lambda$ [2]. The lateral resolution is proportional to the center wavelength and inversely proportional to the numerical aperture (NA) of the optics that focuses the light onto the sample. To image fine details of the sample, a high resolution is needed in both the axial and the lateral direction.
The lateral resolution can be improved by using light with a shorter wavelength or increasing the NA [2]. OCT with a high NA objective lens has also been called optical coherence microscopy [3,4]. With lateral resolutions reaching the micrometer level, sub-cellular structures could be visualized [3]. However, the use of a large NA leads to a limited range where the light is tightly focused, the depth of field (DOF), which is inversely proportional to the square of the NA. The limited DOF is especially problematic for Fourier domain OCT (FD-OCT) where the signal from the full axial range is acquired at once. Thus, most OCT setups use a low NA to capture data that is well-focused over a large depth range. However, these low NA systems cannot reach micrometer-level lateral resolutions.
Fortunately, there are several methods to extend the depth of field. Most hardware-based methods engineer beams with a large DOF, such as Bessel beams [5], or obtain images with different focus depths and combine them after acquisition [6]. Computational DOF extension methods use the overlap between the out-of-focus fields to correct for the defocus [1]. Digital refocusing corrects each en face plane for the defocus by propagating the complex field [1,7]. Interferometric synthetic aperture microscopy (ISAM) uses an inverse scattering model to refocus the whole volume by interpolating in the 3D spatial frequency domain [1,8,9]. These methods can extend the DOF to over an order of magnitude (e.g. 24 times [10]), with the DOF extension being mainly limited by the signal-to-noise ratio (SNR) and the lateral extent of the phase stability.
The lateral resolution can also be improved by reducing the wavelength, with a shorter wavelength leading to a linear reduction of the DOF. As the axial resolution scales quadratically with the center wavelength, reducing the wavelength also improves the axial resolution. Therefore, many high-resolution OCT systems use light sources in the visible (VIS) wavelength range [11,12]. While in the infrared region, the superluminescent diode (SLD) has become the standard light source for spectral domain OCT (SD-OCT), these are not readily available in the visible range. The vast majority of the reported visible wavelength range SD-OCT setups use supercontinuum (SC) lasers [11–13] whose broad spectrum can create a very high axial resolution. However, SC lasers are expensive and the more affordable lasers suffer from high intensity noise. Moreover, the high operating power of the SC laser complicates laser safety requirements and its ultra-broad bandwidth makes the optic design of spectrometers and single mode fiber architecture more complicated.
Superluminescent diodes in the VIS spectrum are available, but the currently available SLDs have a relatively narrow bandwidth. Lichtenegger et al. used a combination of a green (510 nm) and a red (635 nm) SLD together with deep learning and supercontinuum-generated training images on the same system to obtain high-resolution OCT images [14]. However, this still needs a spectrometer to cover the full spectral range, as well as an SC laser to create training data on relevant samples. Khan et al. used an SLD in the blue range (450 nm) for optical coherence microscopy (OCM), but their axial resolution of 12 µm is an order of magnitude above the lateral resolution, strongly limiting the sectioning capability [15]. Thus obtaining both a high axial and lateral resolution is hard with the available VIS SLDs.
Computational improvement of the axial resolution, using deconvolution [16] or spectral estimation methods [17], can be an alternative to using sources with a large bandwidth. Spectral estimation methods can computationally improve the axial resolution in OCT beyond the bandwidth limitation [17,18]. Recently, we presented a fast version of the iterative adaptive approach (IAA) using a recursive scheme and a fast algorithm [18]. Contrary to, for example, the autoregressive (AR) method [17], IAA is parameter free and is able to estimate both the amplitude and the phase of the object. This makes it possible to combine SE-OCT with complex field-based computational methods, such as refocusing [7], inverse scattering algorithms [8] and computational aberration correction (CAO) [9,19].
Here we demonstrate combined high axial and lateral resolution OCT imaging over an extended DOF with a narrow band visible SLD light source. To achieve this, we combine the previously developed IAA-based spectral estimation OCT (SE-OCT) with depth of field extension by ISAM.
First, we apply SE-OCT to improve the axial resolution. The improved axial resolution corresponds to an extrapolation of the interference spectra in the $k$-direction. Instead of using the DFT of the $z$-domain IAA data, the extrapolation in $k$-domain is implemented by using the missing-data IAA (MIAA) [20]. Rather than providing an image in the $z$-domain, MIAA gives the extrapolated data in $k$-domain corresponding to the high axial resolution image. Second, we apply ISAM to extend the DOF. ISAM is a resampling of the OCT spectra in the spatial frequency domain. Third, we apply CAO to correct any remaining aberrations. To test the algorithms and to get insight into the origin of the obtained improvements we have simulated the entire SE-ISAM processing pipeline. The simulations are compared to the experimental data. With MIAA we demonstrate the use of the cost-effective SLDs in the VIS range for high axial resolution OCT while, with the combination with ISAM and CAO, a high resolution is obtained in the lateral direction over a depth range that is much larger than the conventional DOF.
2. Theory
Figure 1 gives a schematic overview of the proposed method, combining spectral estimation and ISAM, applied on the spatial-domain OCT data that is obtained after conventional DFT-based processing. The axial region of interest of the complex-valued OCT image is Fourier transformed into the $k$ direction (Fig. 1(a)) to obtain input spectra for the recursive fast IAA (RFIAA). RFIAA (blue arrow) first normalizes the interference signal and selects the high SNR part (b), after which the actual RFIAA algorithm is applied (Fig. 1(c)). Missing-data IAA (MIAA) (red arrow) not only estimates the spatial domain signal as RFIAA does, but extrapolates the input spectrum (Fig. 1(d)). Then, ISAM is applied (orange arrow) to refocus the image outside the depth of field, resulting in a high resolution in both the lateral and the axial direction. ISAM is implemented via a lateral spatial Fourier transform (Fig. 1(e)), interpolation in $k$-space (Fig. 1(f)), and a 3D DFT (Fig. 1(g)).
In Fig. 1 we propose to first use SE and then ISAM. The opposite order could have an obvious advantage of an improved SNR with ISAM, leading to a higher axial resolution with the application of SE. However, first applying ISAM and then SE has some significant problems on which we will elaborate in section 5.
In the next sections, we first discuss MIAA, including a summary of the previously published RFIAA method. Then we briefly summarize the ISAM theory.
Fig. 1. Schematic overview of the combined SE ISAM processing pipeline. In step (a), the interference spectra are obtained from the low-resolution OCT image. Steps (b-c) show the RFIAA spectral estimation to obtain the complex-valued OCT image with axial resolution improvement. Step (d) added to RFIAA completes MIAA, resulting in the extrapolated interference spectra. Steps (e-g) perform ISAM inverse scattering. The combination of SE and ISAM results in an image with a high resolution in both the axial and the lateral direction. The white dashed lines indicate the edges of the input spectrum for RFIAA spectral estimation OCT. The line plots are illustrations of the spectrum at the center (red-dashed line), with the RFIAA input spectrum in red.
2.1 Spectrum extrapolation with missing-data IAA
2.1.1 Spectral-domain OCT signal
The one-dimensional interference signal in Fourier-domain OCT (in our case SD-OCT) can be described as [21]
where $S_0(k)$ is the source spectral density as a function of wavenumber $k$, $\tilde {a}(z) = a(z)+a^{*}(-z)$ is the combined reflectivity and conjugate reflectivity as a function of propagation distance $z$ in optical path length (OPL). Equation (1) shows that the interference signal, measured as a function of $k$ is the product of the source spectrum and the Fourier transform (FT) of the reflectivity $\tilde {a}(z)$. The reflectivity $\tilde {a}(z)$ is usually estimated with an inverse FT resulting in an OCT reflectivity that is the convolution of the reflectivity with the inverse FT of the source spectrum, the latter of which acts as the axial PSF.In OCT, the interference signal is measured on a discrete grid and can be approximated as [18]
where $\mathbf {y}_g$ is an $N_g\times 1$ vector containing the given normalized interference signal $I(k_n)/S(k_n)$ at discrete wavenumber $k_n$ and $a(z_l)$ is the discretized reflectivity. The vector2.1.2 Iterative adaptive approach
Obtaining the reflectivity $a(z_l)$ from the measured spectrum is equivalent to a spectral estimation problem [17,18], which can be solved with a variety of methods. Recently, we showed that the non-parametric iterative adaptive approach (IAA) [22] can significantly improve the axial resolution with respect to the conventional discrete Fourier transform (DFT) reconstruction [18].
In brief, IAA estimates the reflectivity $a(z_l)$ with a weighted least squares solution of Eq. (2) as
The solution of Eq. (4) can be written as
The sample reflectivity $a(z_l)$ is initialized as the DFT of $\mathbf {y}_g$ with zero-padding, which is equivalent to initializing the data covariance matrix with the identity matrix, $\mathbf {R}_g=\mathbf {I}$. By iterating between Eqs. (5) and (6), the estimate for the (complex) reflectivity $a(z_l)$ is refined. Usually, 10 iterations are sufficient for convergence [22].
2.1.3 Missing-data IAA
IAA estimates the spatial-domain reflectivity $a(z_l)$, however, for the application of ISAM, the spectral data corresponding to the high axial resolution image needs to be estimated. This is an extrapolation of the interference spectrum outside the range where it is measured, i.e. this data is missing from the measured data.
When missing data problems occur, the IAA algorithm can be used for missing data recovery following a two step procedure as described in [20]. MIAA can cope with arbitrary missing data patterns, for uniform or nonuniform sampling, interpolation as well as extrapolation of data sequences. In our case, we consider the interference spectrum at wavenumbers $k$ where it is measured (to be more precise, the high SNR part of the measured interference spectrum) as the given data. For the application of ISAM, we aim to estimate the interference spectrum at wavenumbers outside the given spectral range (left-hand side and right-hand side data extrapolation). This data is not physically measured and will be indicated as missing data. First, IAA is applied on the given data set for the estimation of the $z$-domain reflectivity parameters and the computation of the relevant covariance matrices. Then, extrapolated data (missing data) are estimated using a linear minimum mean squared error (MMSE) estimator.
The MMSE MIAA data extrapolation considers the data sequence of the interference spectrum of the form
In the first step of MIAA, the reflectivity $a(z_l)$, $l=0,1,\ldots, L-1$ has been estimated using IAA on the given data, iterating Eqs. (5) and (6). At the missing data recovery step, a general linear estimator of ${\bf y}_m$ is considered, [20],
The matrix ${\bf T}$ that gives MMSE estimation of ${\bf y}_m$ is given by where ${\bf R}_{mg}$ is the $N_m\times N_g$ cross-covariance matrix between the missing or extrapolated data and the given data. ${\bf R}_{mg}$ can be estimated using known quantities in a similar way to the autocovariance matrix in IAA as where ${\bf f}_m(z_l)$ is the $N_m \times 1$ vector with Fourier components at $z_l$ corresponding to the missing wavenumbers $k_n$ equivalent to Eq. (3). Combining Eqs. (9) to (11), we obtain the estimation of the missing data2.1.4 Computational implementation of MIAA
Here, we consider the application of MIAA on OCT spectral data that is uniformly sampled in the $k$-domain. Moreover, the reflectivity coefficients $a(z_l)$ in Eq. (5) are estimated on a uniformly sampled $z$-domain grid. As a result, the components of the Fourier vector in Eq. (3) are of the form
with $l$ and $n$ integers that go from 0 to $L-1$ and $N_g-1$ respectively, implying a discrete space Fourier representation for the Fourier vector.The specific structure of the given data as well as the Fourier vector, Eqs. (7) and (8) respectively, are such that the given data consists of a continuous data segment ${\bf y}_g$. This allows for the use of the fast IAA (FIAA) method for the estimation of the reflectivity coefficients $a(z_l)$. Thus Eqs. (5) and (6) are actually implemented in a computationally efficient way using fast Toeplitz matrix algebra and the Fast Fourier Transform (FFT) [18,23]. The computational complexity of a single iteration of the FIAA algorithm is given by
which compares favorably against the ${\mathcal {O}}\left ( N_g^3 + N_g^2L \right )$ complexity required by the brute force approach.Moreover, the missing data consist of continuous data segments ${\bf y}_{mL}$ and ${\bf y}_{mR}$ which are actually on the left-hand side and the right-hand side of the interpolated spectral data, respectively. Thus, we implement Eq. (12) using fast Toeplitz matrix operations and the FFT, resulting in an additional computational cost of
Consequently, the overall computational cost of the proposed fast MIAA for the recovery of the missing OCT data is given byB-scan OCT imaging is performed by processing consecutive A-scans as columns in an image matrix. Although these columns can be processed independently to produce the corresponding sequence of depth profiles, a warm start initialization procedure can be applied that drastically reduces the required amount of iterations. We use the data covariance matrix of the previously processed A-scan for the initialization of the currently processed A-scan as it is expected that successive A-lines have only a slight variation to each other as they partially probe the same sample structure. We called this method the recursive fast IAA (RFIAA) [18,24] and this approach generates results without any substantial loss in performance using only 2 iterations. When RFIAA is used in place of FIAA, the computational complexity of the proposed fast recursive MIAA for the recovery of the missing OCT data is given on the average (per A-scan) by
2.2 Interferometric synthetic aperture microscopy
The 3D measured signal in point-scanning OCT can be described as the convolution of a space-variant complex PSF $h(x,y,z;k)$ with the scattering potential $f(x,y,z)$ of the sample [1] as
where $x$ and $y$ are the lateral coordinates of the beam location, $z=0$ is the focus depth, and $k$ is the measured wavenumber corrected for the refractive index of the sample.By taking the lateral Fourier transform and using asymptotic approximations for near-focus and far-from-focus cases, Eq. (18) can be rewritten as [1,25]
where $H(k_x,k_y;k)$ is the space-invariant optical transfer function, and the integral is the Fourier transform in the $z$-direction of $f(k_x,k_y,z)$ with wavenumber For Gaussian beams $H(k_x,k_y;k)$ is relatively smooth and acts as an amplitude optical transfer function. Thus, the scattering potential in $k$-space can be approximated asInterpreting ISAM in the $k$-space description, the OCT signal for a wavenumber $k$ is obtained along the Ewald sphere in the $(k_x,k_y,k_z)$-space with its center in the origin and a radius of $2k$, the factor 2 accounting for the backscattering geometry [26]. Rather than assuming that $k_z=2k$, which is done when the axial reconstruction is considered independent from the lateral spatial frequency, ISAM places the data at its true $k_z$ coordinate as given by Eq. (20). After interpolation to a linear grid in $k_z$, see Fig. 1(f), the refocused image with depth-invariant resolution can be obtained by taking a 3D inverse DFT of $f(k_x,k_y,k_z)$, see Fig. 1(g).
3. Methods
3.1 Experimental setup
Figure 2 shows a schematic overview of the custom build high-resolution OCT setup. Light from the fiber-coupled green superluminescent diode (EXS210118-01, Exalos) is coupled into a 50:50 wideband fiber coupler (TW560R5A2, Thorlabs), that distributes the light to the reference and sample arm. The reference arm consists of a collimator lens (AC256-050-A, Thorlabs), a mirror, and an iris to control the reference arm light power. A fiber-based polarization controller (FPC560, Thorlabs) in the reference arm is used to align the reference arm polarization with that from the sample arm, thus optimizing the interference signal. In the sample arm, a collimator lens (AC080-020-A, Thorlabs) gives a Gaussian beam with a waist of 3.28 mm. Two galvo mirrors (RTA-AR180, Newson, Belgium) are placed around the back focal plane of a scan lens (CLS-SL, Thorlabs) that is followed by a matched tube lens (ITL200, Thorlabs) for telecentric scanning and an objective lens (10x Plan Apochromat, Mitutoyo) with an NA of 0.28 and a working distance of 34.0 mm. The 11.2 mm diameter aperture of the objective lens is almost completely filled with the beam (expanded by the scan lens and tube lens) that has a waist diameter of 9.4 mm. The ratio between the scan angle and lateral displacement of the focus was experimentally calibrated using a resolution test target (R1DS1N, Thorlabs) as a sample. The obtained values give a lateral FOV of 1.28 mm $\times$ 1.35 mm for the $\pm 5 ^{\circ }$ scan range of the galvo mirrors, which is close to the theoretical FOV of 1.22$\times$1.22 mm$^2$. The light from the reference and sample arm is recombined with the fiber coupler and guided to the spectrometer.
In the custom-build spectrometer, the light from the fiber is collimated with a collimator lens (F220APC-532, Thorlabs), expanded by a beam expander (ACN-127-A and AC508-300-A), and projected on a 50.8 mm diameter volume phase holographic grating with 1800 lines/mm (Wasatch Photonics, USA). The dispersed beam is imaged on a 6144-pixel line-scan camera (raL6144-80km, Basler, Germany) with a focusing lens consisting of two identical achromatic doublets (AC508-750-A, Thorlabs). The large beam diameter on the grating combined with the optics design based on simulations with Zemax ensured a good spectral resolution.
Fig. 2. Schematic drawing of the experimental setup. SLD: superluminescent diode, BS: fiber beamsplitter, PC: polarization control, CL: collimator lens, Ir: Iris, M: mirror, GM: galvo scan mirrors, SL: scan lens, TL: tube lens, OL: objective lens, SH: sample holder, BE: beam expander, G: grating, FL: focusing lens, CA: Line-scan camera.
The chirp values for $k$-linearization were obtained from the difference in the unwrapped phase from two averaged A-scans with a mirror on either side of the zero-delay in the sample arm [27]. The dispersion mismatch was determined from a measurement of a single mirror with the phase fitted with an 8-order polynomial. The fitted phase deviation was used for dispersion compensation by multiplying the interference signal with a complex exponential of the phase difference [28]. To obtain sufficient SNR over the entire spectral width of the spectrometer, the $k$-linearization and dispersion mismatch calibration was performed using light from a fiber-coupled supercontinuum light source (NKT EUL-10, NKT photonics). The depth sampling density was determined from a linear fit on the peak locations of 32 averaged A-scans with a mirror position translated over a total range of 4.4 mm. The maximum imaging depth of 8.27 mm in air was then used to determine the linearized $k$-sampling density at 190 m$^{-1}$. A diode laser at 532 nm (CPS532, Thorlabs) was used as a reference wavelength to obtain the physical $k$-values for the linearized spectrum. The full spectrometer range is 48.7 nm (from 488.3 nm to 537.0 nm) with a sampling density of 7.9 pm per pixel. Using a fit on the sensitivity decay [29], the effective spectral resolution was determined to be 13.4 pm. The sensitivity decay is 6 dB at 4.1 mm ($z_{max}/2$) and 13 dB at 6.2 mm ($3z_{max}/4$). For the OCT measurements with the SLD, only 3072 of the 6144 pixels were used between 501.3 nm and 525.6 nm as the intensity is insignificant outside this region. The measured FWHM of the light source intensity spectrum is 6.5 nm. The FWHM of the axial PSF in air was measured to be 11 µm.
The camera was triggered per B-scan on the framegrabber (PCIe-1433, National Instruments) by the galvo mirrors, the angular sweep speed being adapted to the line rate of the camera. The data acquisition was done with a custom script in Python 3.7, operated with Anaconda Spyder on a desktop computer. Also, all simulations and pre-processing of the data was done in Python 3.7.
3.1.1 Samples and acquisition settings
Two samples were used to evaluate the performance of the proposed method. The first test sample consisted of TiO$_2$ powder (Sigma Aldrich) in gelatin (Dr. Oetker, the Netherlands). A droplet of a diluted suspension of TiO$_2$ particles in water was added to the heated water-gelatin mixture. The mixture was poured into a custom mount, covered by a coverslip to create a flat top surface, and cooled down to room temperature. The sub-resolution particles were used to characterize the 3D resolution.
The second sample was a leaf disc that was punched out of a lettuce leaf. The leaf disc was water infiltrated by putting it in a syringe with water and lowering the pressure. When the pressure is released, the gas-exchange cavities are filled with water, reducing the refractive index contrast [30]. The leaf disc was mounted in water between two coverslips, at a sufficient distance from the coverslip to avoid saturation artifacts from the coverslip reflection.
From each sample a volume of $512\times 512$ scan lines was obtained over a lateral area of 0.225$\times$0.225 mm$^2$, giving a lateral sampling of 0.44 µm. The exposure time was set to 80 µs, at a line acquisition rate of 11.8 kHz. The large exposure time was needed because the power on the sample was measured to be only 170 µW, probably due to a limited coupling efficiency between the source and the fiber coupler/splitter.
3.2 OCT signal simulations
Three-dimensional OCT data was simulated by combining 1D OCT spectrum simulations [18,21] with the intensity and phase of a Gaussian beam. The waist in the focus of the Gaussian beam was chosen based on the experimental beam width in the back-focal plane of the objective and the propagation to the relevant depth. The aperture edge is chosen where the Gaussian beam intensity is $1/\text {e}^2$, giving an effective NA of 0.235 with the 9.4 mm waist diameter. Then the waist in focus can be calculated as
where $\lambda _c$ is the center wavelength. The OCT interference signal for a discrete set of reflectors isThe OCT interference spectrum for each lateral position $(x,y)$ is then calculated as
where $E_r(k)=\sqrt {S_0(k)}$ is the reference beam field. The simulated OCT image is obtained by taking the inverse DFT of $I(x,y,k)$ from $k$ to $z$. Subsequently, complex Gaussian white noise is added with a standard deviation that depends on the intensity of the scatterers and the desired SNR.3.3 Data processing
3.3.1 Conventional OCT processing and phase correction
The OCT spectral data was first processed in the conventional way by subtracting the reference spectrum, interpolating to a linear grid in $k$-domain, multiplicating with a dispersion correction vector, and performing an inverse DFT. The reference spectra were obtained by laterally averaging the interference signal of a B-scan, excluding the spectra with saturation.
As ISAM requires lateral phase stability, the experimental data was corrected for phase drift during the measurement using a coverslip interface as a reference. The phase of the coverslip was determined by putting the image outside the signal of the coverslip at the positive depth, everywhere to 0, taking the DFT of the complex-valued OCT data in the axial direction, unwrapping the phase, and applying a linear fit to the unwrapped phase. The difference between the fit and a reference slope was corrected for by multiplying the entire spectrum with a complex exponential, such that after correction, the coverslip became a flat, horizontal plane with a constant phase.
To reduce memory requirements and computational complexity, the complex-valued image over an axial region of interest of 200 pixels (0.81 mm in water) was selected for further processing.
3.3.2 RFIAA and MIAA processing
The axial DFT of the phase-corrected image was used as input for the RFIAA spectral estimation algorithm. The interference signal was normalized using the reference spectrum. For the experimental data, this was obtained from the lateral average of the absolute value of the interference signal of the volume. Subsequently, RFIAA was applied on a 128-pixel-long central part of the spectrum with the highest intensity. The RFIAA reconstruction grid length $L$ was 800 pixels, four times the original 200-pixel grid length, giving an increase in spatial sampling by a factor of 4 with respect to the original image. The number of iterations for RFIAA was set to 10 for the first line and 2 for subsequent A-lines that were initialized with the covariance matrix of the previous A-line. For MIAA, the RFIAA reconstruction grid length $L$ was chosen eight times longer than the original grid length, such that number of samples $N_m+N_g$ is four times the length of the 200-pixel long original image. After RFIAA, the MAP estimation was applied, followed by an FFT, a circular shift, and the selection of $N_m$ values to obtain the extrapolated spectra.
To reduce side lobes in the axial PSF, the edges of the MIAA extrapolated spectrum were tapered with a 200-pixel-long squared cosine window on both sides, bringing the intensity smoothly to 0. After the lateral Fourier transform, the high lateral frequencies, mainly containing side lobes and noise, were apodized with a circular window with a squared cosine-shaped edge with an inner radius of 110 pixels (3.1 µm$^{-1}$) and an outer radius of 210 pixels (5.9 µm$^{-1}$). While the side-lobes were significantly suppressed, the effect on the main lobe width for the experimental data was barely noticeable, while the effect on the simulation data was moderate.
RFIAA and MIAA were implemented in MATLAB (R2020a). The processing was done on a Dell Precision 5820 desktop computer with an Intel Xeon W-2223 CPU and 32 GB RAM. The processing time for the volume with $512\times 512$ A-lines was around 113 s for RFIAA and 181 s for MIAA, 0.22 s and 0.35 s per B-scan respectively.
3.3.3 ISAM data processing
ISAM requires the focal plane depth and the sample refractive index as input parameters. The focal plane depth index was obtained by visually inspecting where the image appeared most sharp. The used refractive index was set to $n=1.33$, the refractive index of water at 900 nm wavelength and room temperature [31]. The sampling period of $k_z$ was twice that of $k$, such that $k_x=k_y=0$, $k_z=2k$ coincides with the original grid. The grid of $k_z$ was extended such that the $k$-space contained all the interpolated data. After MIAA and ISAM the $z$ grid had a length of 1029 pixels in depth, resulting in an axial sampling density of 0.79 µm per pixel.
ISAM processing was implemented in MATLAB and took around 60 s per volume of $512\times 512$ A-lines. The processing time for the conventional DFT reconstruction with ISAM was the same, as zero-padded data was used to obtain the same axial sampling density as with MIAA.
3.3.4 OCT resolution analysis
The 3D resolution was quantified from the images of the simulated and measured point scatterer objects. The positions of the scatterers were automatically determined using a 2D local maximum detection algorithm [32] in the $xy$-plane and the $xz$-plane. For the local maximum detection, the image is first divided by the square root of the laterally averaged intensity, that is axially smoothed with a 40-pixel standard deviation Gaussian kernel. This reduces the intensity difference between the in-focus and out-of-focus regions and enables the use of a global threshold in the final scatterer selection. Then, the difference between two local average filters with a small (4-pixel diameter) and a large (8-pixel diameter) kernel is calculated. The local maximum position is obtained where the value of this difference is the maximum over a 12-pixel-diameter kernel region. A global threshold of 1/2000 of the maximum intensity in the image excludes noise from the detected local maxima. The 3D location of the scatterers is where the 2D local maxima of the orthogonal planes overlap.
For each 3D scatterer position a $21\times 21\times 21$ pixel ($16.6\times 9.2\times 9.2$ µm$^{3}$) subvolume is selected. If the maximum intensity inside the volume lies outside the $3\times 3\times 3$ center pixels, the volume is excluded as it contains another scatterer with a higher intensity that would contaminate the fit. For DFT and RFIAA without ISAM, the locations from the ISAM images were used with a $31\times 201\times 201$ pixels ($24\times 88\times 88$ µm$^{3}$) subvolume, to capture the full defocused PSF.
To quantify the resolution, a 3D Gaussian function was fitted to the intensity distribution in the selected subvolumes using the least square curve fitting function lsqcurvefit. The standard deviation in every direction was then used to calculate the full width at half maximum (FWHM) which serves as a resolution measure. The SNR for each scatter was determined by dividing the peak intensity of the Gaussian fit by the variance of the noise amplitude in a manually selected area of the image without signal.
3.3.5 Computational adaptive optics
Computational adaptive optics (CAO) was implemented using the sub-aperture correlation method [33] implemented on en face planes of the OCT data. The complex-valued OCT en face images were Fourier transformed in the lateral direction to obtain the field in the pupil plane. The field in the pupil plane was shifted such that the lateral center of mass was in the center [34], after which it was transformed back to the spatial domain. These en face images were divided into four square areas, each of which was independently corrected to also address shift-variant aberrations. The image from each area was lateral Fourier transformed to obtain the field in the pupil plane. In the pupil plane field, 45 circular sub-apertures were defined on a $7\times 7$ grid (with the four corners excluded) within a circular aperture with a radius of 1 µm$^{-1}$. With sub-aperture radii equal to the lateral spacing between the sub-aperture centers, the sub-apertures were partially overlapping. After inverse Fourier transformation, the magnitude images of each sub-aperture were cross-correlated with the magnitude image of two randomly chosen sub-apertures. The obtained lateral shift between the images is proportional to the difference in the local slope of the aberration between both sub-aperture positions. The obtained slope differences were fitted with the $x$ and $y$ gradient of 12 Zernike polynomials (from second to fourth order) [33]. By using two randomly chosen sub-apertures instead of the center aperture as reference, the sensitivity to speckle was reduced [34,35]. Magnitude OCT images were used rather than intensity OCT images to increase the sensitivity to the entire object structure rather than a few high-intensity peaks.
As the speckle pattern from non-overlapping sub-apertures is independent, images with strong speckle structures can give erroneous shift estimates, even when using the more robust correlation with multiple randomly chosen apertures [35]. As the aberrations are expected to change smoothly over depth, the robustness was further increased by using a weighted moving average in the depth direction over the Zernike coefficients. A 10-pixel sigma Gaussian kernel, combined with the inverse of the least squares fit error for weighting, ensured both smoothness and a lower weight for unreliable coefficients. The CAO was completed by phase conjugation in the pupil plane using multiplication by a complex exponential to obtain a flat phase profile. The aberration-corrected en face OCT image was then obtained by inverse Fourier transformation. As ISAM removes most of the defocus, the magnitude of the aberrations was limited and it was not needed to iterate the procedure.
CAO was implemented in Python 3.7 and took around 4 s per en face image (1 s per sub-image).
4. Results
4.1 Point scatterer sample OCT imaging
Figure 3 shows the OCT images of the gelatin phantom filled with TiO$_2$ particles. The TiO$_2$ particles are well visible in the OCT images as they have a high intensity. The B-scans show a band of high-intensity background near the focal plane. They come from impurities in the gelatin that appear in the focal plane depth range because of the high light intensity.
The DFT method shows vertical stripes near the focus in the $xz$-plane, Fig. 3(a) due to the mismatch between the axial and the lateral resolution. The en face image in the focal plane, Fig. 3(b), shows many spots at a resolution around the diffraction limit. Outside the small focal region, Fig. 3(c), large blurred spots are visible in the $xz$-plane and the out-of-focus en face image.
Fig. 3. OCT images of the TiO$_2$ sample processed with the different methods: (a) shows B-scans, (b) shows sections of the en face image in focus (blue line in (a)), and (c) shows the en face image 82 µm above focus (red line in (a)). The white arrows indicate the scatterers that are shown in Fig. 4(a-f) in more detail. A 3D rendering of the data in the figure is shown in Visualization 1.
RFIAA enhances the axial resolution, giving almost isotropic-shaped spots in the focal plane. The in-focus en face image (b) shows fewer scatterers as the optical depth sectioning due to the spectral estimation is improved. However, as RFIAA does not account for the defocus, there is still a large lateral blur outside the focal range as seen in Fig. 3(a) and (c).
With DFT+ISAM, the signal outside the focal region is effectively refocused, giving narrow vertical stripes (Fig. 3(a)) or points in the en face image (Fig. 3(c)). The out-of-focus lateral resolution is similar to that in focus, as expected from the correct implementation of ISAM. Due to the poor axial resolution, the en face image shows many scatterers and the B-scan image has limited quality.
Combining MIAA and ISAM both improves the axial resolution and extents the depth of field. The axial resolution and optical sectioning in focus are similar to that of RFIAA. Outside the focal depth range, ISAM refocuses the signal, which is visible in the narrow spots in the $xz$-plane Fig. 3(a) and the small dots in Fig. 3(c). In the en face images with MIAA+ISAM, fewer scatterers are visible than with DFT+ISAM due to the improved axial resolution. Some background signal is visible around the scatterers at the location of the original defocused spot, which we attribute to residual phase noise in the MIAA reconstruction at low SNR, giving incoherent signals that cannot be refocused with ISAM. As the intensity is much lower than that of the scatterer, the impact on the image quality is limited.
The imaging results in Fig. 3 qualitatively show that MIAA+ISAM successfully combines the axial resolution improvement of IAA-based SE-OCT with DOF extension using ISAM. The next section analyzes these results quantitatively.
4.2 OCT resolution analysis
Figure 4 shows the results from the resolution analysis of the OCT images of the TiO$_2$ sample. Figure 4(a-f) shows two exemplary point scatterers reconstructed with the four different methods and fitted with a 3D Gaussian function. The in-focus scatterer (Fig. 4(a-c)) is imaged with relatively high resolution for all four methods and the line plots indicate that the PSF is well described by the Gaussian model. The axial resolution improvement with MIAA+ISAM is, with an axial resolution of 1.45 µm, a factor 5 better than without MIAA. Outside focus (Fig. 4(d-f)) DFT and RFIAA without ISAM give very blurry spots because of interference between the signals of several defocused scatterers and because the PSF width is larger than the chosen small volume. This also gives unreliable Gaussian fitting, see Fig. 4(e-f), thus we excluded these methods from the further data analysis. MIAA+ISAM is, with a resolution of 4.0 µm, a factor 1.8 better than without MIAA. The smaller improvement in axial resolution compared to in focus is caused by the lower SNR in the out-of-focus region.
Fig. 4. The lateral and axial OCT resolution for the various computational methods. The top and middle row show an example point scatterer at focus (white arrow in Fig. 3(b)) and 82 µm above focus (white arrow in Fig. 3(c)). (a,d) show the $xz$ cross sections, (b,e) the lateral and (c,f) axial line cross sections at the white dotted line. The circles indicate the OCT intensity and the solid line is the corresponding 3D Gaussian fit. (g) shows the SNR as a function of depth from focus, the experimental data being grouped per 20 µm depth range. (h-i) show the resolution in the axial and the lateral direction for experimental and simulated data. In (i) the lines for the simulated data without (DFT and RFIAA) and with ISAM (DFT+ISAM and MIAA+ISAM) overlap.
The 809 automatically detected points that were at the same position ($\pm 1$ pixel in the axial direction) for both DFT+ISAM and MIAA+ISAM were used for the analysis. Simulations of 48 point scatterers divided over 4 vertical columns that span a 450 µm depth range were used to further interpret the experimental results. Here, also DFT and RFIAA without ISAM were included as the simulated scatterers were well-separated.
Figure 4(g) shows how, for all methods, the SNR peaks in focus and decreases further away from focus. For the experimental DFT/MIAA + ISAM data, the SNR in focus is on average around 60 dB with peaks to 80 dB, from where the SNR decreases to 40 dB at 150 µm from focus. The error bars indicate the lowest and highest SNR within each 20 µm depth interval. The large variation, especially in focus can be attributed to variation in the scattering amplitude of the TiO$_2$ particles. Both DFT and MIAA with ISAM are in good agreement with the simulations and have a high SNR outside of focus due to the ISAM refocusing. MIAA+ISAM has a slightly lower SNR at the edges. The simulations for DFT and RFIAA show that in focus, the SNR is similar, outside focus it drops strongly because the signal is smeared out over a large lateral range and not refocused with ISAM. The SNR with RFIAA is a bit lower due to a bit higher noise variance. The SNR of the DFT simulation drops to 10 dB at 200 µm from focus. The coherence of the low-intensity signals remains however sufficient for applying ISAM, which then improves the SNR outside focus with 15 to 30 dB.
Figure 4(h) shows the axial resolution as a function of the depth from focus. As expected, DFT+ISAM gives a depth-independent axial resolution of around 8 µm, slightly below the bandwidth-limited resolution of 8.1 µm. Simulations show a similar depth-independent axial resolution. ISAM improves the axial resolution with about 1 µm because it maps high lateral frequencies to lower axial wavenumbers $k_z$, giving an effective broadening of the spectrum (see Fig. 1(g)). This effect is especially significant for sources with a narrow bandwidth as the relative curvature in $k$ space is large.
MIAA+ISAM gives a significant improvement in the axial resolution to around 1.4 to 2 µm in focus and increasing to 4.5 µm at 150 µm distance from the focus. The variation in axial resolution with depth can be explained by the variation of the SNR with depth. An SNR before MIAA above 50 dB gives a factor 5 improvement, an SNR of 30 dB gives a factor 3 improvement and an SNR of 10 dB reduces the improvement to a factor 1.5. The improvement in axial resolution, especially for SNR above 30 dB, brings it closer to being similar to the lateral resolution. The simulations show that for RFIAA and MIAA the in-focus axial resolution of 0.95 µm is slightly better than the experimental resolution. The poorer experimental axial resolution could be attributed to the background scattering of the gelatin, which reduces the local SNR. The 0.95 µm is close to the theoretical resolution of 0.85 µm, assuming a uniform SE bandwidth extension. Outside focus, RFIAA without ISAM gives a slightly worse axial resolution than MIAA+ISAM, because of the spectrum broadening with ISAM.
Figure 4(i) shows the OCT lateral resolution analysis. Both DFT and MIAA in combination with ISAM give a depth-independent resolution of around the diffraction limited value of 0.85 µm. This shows that MIAA sufficiently preserves coherent phase information for applying ISAM. The variation in resolution can be caused by a variation in scatterer size or by wavefront aberrations. These aberrations, which are also the cause of the large deviation at depths below focus, can be computationally corrected [9,35], as we will show in section 4.4. As the deviation is similar for DFT and MIAA, it is not caused by MIAA. Simulations show that without ISAM, the nominal DOF is only about 25 µm. DFT and MIAA in combination with ISAM give a strong improvement in the resolution outside the focal region. Their depth-independent lateral resolution is around 0.74 µm over the full simulated depth range of 450 µm. The simulated resolution is at the lower limit of the experimental results because the simulations do not include the physical aperture of the objective lens. However, due to the lateral apodization during processing, the difference is limited.
The results on the lateral resolution clearly show the effectiveness of ISAM to extend the depth of field, also after spectral estimation. The shown extended DOF is around 400 µm (a factor 16 extension) for experimental data, and at least 450 µm (a factor 18 extension) for the simulation data. Note that the DOF with ISAM is more constrained by the SNR limiting the axial resolution than the remaining lateral blurring.
4.3 Plant leaf computational OCT imaging
Figure 5 shows the results for a more realistic object, a lettuce leaf. The comparison between DFT and MIAA, both in combination with ISAM clearly shows the improvement in axial resolution and image quality that our method achieves. While the plant cells in Fig. 5(a) barely can be distinguished with DFT+ISAM, they are clearly visible with MIAA+ISAM. The improved optical sectioning ability is also clear from the en face images, Fig. 5(c), where MIAA+ISAM better visualizes the open leaf structure. This allows for clearer visualization of the plant cells and sub-cellular structures. An example of the improved visualization of open structure with MIAA+ISAM is shown in the zoom-in. The object has a better contrast between the water-filled intra-cellular space and the cell walls. The cell walls also appear narrower with MIAA+ISAM, probably because a tilted cell wall is laterally smeared out due to the poor axial resolution of DFT+ISAM.
Fig. 5. OCT imaging of the lettuce leaf. (a) An $xz$-cross section of the 3D leaf reconstruction, with (b) the zoom-in at the yellow square for all four methods. (c) The en face image at the white dashed line in (a), the white dashed line indicates the intersection between (a) and (c). The bottom images show the zoom-in at the blue square. (d) The en face image at the bottom of the epidermal cell layer indicated by the green dashed line in (a), for all four methods. The green dashed line indicates the intersection between the images of (a) and (d). (e) The A-scan at the vertical white line in (a) for all four methods, the points indicate OCT intensity, while the solid line is a 2-peak Gaussian fit through the intensities. A 3D rendering of the data in this figure is included in Visualization 2.
Figure 5(b) shows a close-up of the OCT B-scan (a) for all four methods. The large asymmetry between axial and lateral resolution makes it hard to distinguish the cell walls in both DFT-based methods. RFIAA already makes the structure much clearer. However, the lateral blurring and intensity fluctuations caused by the interference of blurred signal significantly decrease the image quality, especially at the top half. By combining the axial resolution improvement of spectral estimation with the out-of-focus lateral resolution enhancement of ISAM, the cellular structures are also clearly imaged in the top half of the inset.
The effect of ISAM is even clearer from the en face images in Fig. 5(d), which are located around the bottom of the epidermal cell layer, the green dashed line in Fig. 5(a). Without ISAM, the image is significantly blurred and does not reveal any plant structure. ISAM refocuses that into a clear image of the plant tissue structure. However, with DFT+ISAM the structures of different depth layers are merged. For example, the structures indicated with the green arrows have their peak intensity in a plane 6 µm higher up, the structure at the red arrow shows up in the MIAA+ISAM image 2.5 µm up and the structure indicated with the white arrow is located 4 µm lower. MIAA+ISAM successfully sections out the signal from a narrow depth region and refocuses it with a good lateral resolution.
Figure 5(e) shows the A-line profiles at the vertical white line in Fig. 5(a). With DFT+ISAM in Fig. 5(a) the shape of the cell is not clear because of the poor axial resolution, while MIAA+ISAM gives a clear image. The solid lines in Fig. 5(e) correspond to a 2-peak Gaussian fit of the A-line segments. MIAA+ISAM gives an FWHM of the peaks of 1.8 µm and 2.1 µm respectively, compared to 7.2 µm and 6.4 µm for DFT+ISAM, a factor 3 to 4 improvement. Without ISAM, the top layer is not clearly visible, probably because of destructive interference of defocused signal. The small shift in peak location with respect to the ISAM methods could be the result of focusing the curved structure.
The lettuce leaf OCT imaging clearly shows that MIAA+ISAM does not only work for ideal point scatterers, but also for relevant biological samples with a medium level of sparsity. This proposed method gives a clear improvement in axial resolution and image quality, resulting in good-quality images with the cost-effective green SLD light source that has a bandwidth (FWHM) of only 6.5 nm.
4.4 Computational adaptive optics OCT
Though ISAM refocuses the signal outside the focal plane, there can be other aberrations and residual defocus. As MIAA+ISAM preserves the phase of the (complex) signal, computational adaptive optics (CAO) can be applied after MIAA+ISAM reconstruction. Figure 6(a-c) shows the results with CAO for the sample with the TiO$_2$ particles. The most significant Zernike coefficients, for defocus $Z_{0}^{2}$ and primary spherical aberration $Z_{0}^{4}$, are displayed in Fig. 6(a) as a function of depth. Most of the other Zernike coefficients varied around 0. Especially the defocus aberration is significant with a coefficient up to 2 radians below the focus. Note that the residual defocus is not due to wrong ISAM implementation as the aberration above focus is close to 0 for a large depth range. In absence of speckle, the error in the estimated coefficients is limited. Still, a few outliers of the estimations (blue dots) are visible, which the weighed moving average (red line) removes.
The en face plane 190 µm below focus (Fig. 6(b)) clearly shows the improved resolution or sharpness of the image after CAO. The improvement is quantified using 3D Gaussian fits on 983 and 1075 particles, respectively. The histogram of the lateral FWHM (Fig. 6(c)) shows that the aberrated spots with an FWHM above 1.2 µm are corrected, while the peak value remains at the diffraction-limited resolution. This correction improves the mean lateral FWHM with 8% to 0.87 µm, but the improvement for the aberrated spots is of course much more significant. Most improvement is obtained for depths below focus, pushing the large spot sizes that were seen in Fig. 4(i) back to diffraction-limited values.
Fig. 6. OCT imaging results for computational adaptive optics (CAO) applied after MIAA+ISAM processing. (a) and (d) show two Zernike coefficients for the top right corner of the FOV as a function of depth, the blue dots being the estimates per plane and the red line the weighed moving average. (b) shows an en face image 190 µm below the focus of the TiO$_2$ particles in gelatin sample before and after CAO. (c) shows the histogram of the lateral resolution before and after CAO. (e) shows an en face image of the leaf at 24 µm below focus before and after CAO, (f) the inset from (e), and (g) the line-plot along the green line in (e). Visualization 3 shows more CAO en face images.
The lettuce leaf sample has more severe aberrations than the sample with the TiO$_2$ particles. The large defocus aberration, shown in Fig. 6(d), could potentially be reduced by optimizing ISAM settings, but that effect would never fully compensate for the defocus as the aberration is not fully linear as a function of depth. Moreover, manual optimization of the focus depth for ISAM did not give a better image over the full depth range. For the plant leaf, the variation in estimated Zernike coefficients from plane to plane is much larger than with the TiO$_2$ sample due to the presence of speckle (note the different scale in Fig. 6(a)). The weighed moving average (red line) reduces the variations and increases the robustness of the CAO. At depths of >100 µm below focus, the image quality was too poor to have reliable aberration estimates, and even with manually optimized defocus the image quality remained poor. This may be caused by the accumulation of phase errors in the heterogeneous medium.
CAO clearly improves the sharpness and quality of the en face image (Fig. 6(e)). The cell walls are imaged much sharper and subtle structures are much clearer visible, as shown in the inset of Fig. 6(f). The line plot in Fig. 6(g) shows three peaks corresponding to cell walls, whose FWHMs improve with a factor 2 from 3.5 µm to 1.6 µm on average.
The results in Fig. 6 show that the data after MIAA+ISAM is of sufficient quality to apply CAO for further optimization of the resolution.
5. Discussion
In this work, we combined axial resolution improvement using SE-OCT with ISAM for depth of field extension.
In this discussion, we first consider our motivation for our data processing by first applying MIAA followed by ISAM. Next, we discuss some limitations of our approach and explore future opportunities.
Instead of applying spectral estimation on the original spectra (in the $k$-direction), applying it to the ISAM-processed data (in $k_z$-direction) could have an important advantage. ISAM improves the SNR of the out-of-focus signal, which would then result in a better axial resolution with SE-OCT. There are several ways in which ISAM and SE can be combined. We discuss here two combinations:
- 1. First applying ISAM, then applying SE from the $(x,y,k_z)$-domain to the $(x,y,z)$-domain.
- 2. First applying ISAM interpolation, leaving the interpolated data in the $(k_x,k_y,k_z)$-domain, then applying SE from $(k_x,k_y,k_z)$-domain to the $(k_x,k_y,z)$-domain, followed by a lateral inverse DFT to the $(x,y,z)$-domain.
The first option would be the most obvious way, but there is a significant complication with performing the image enhancement in this order. As ISAM maps the signal at higher lateral spatial frequencies $k_x, k_y$ to lower $k_z$, the spectrum in $k_z$ will be broader than the spectrum in $k$ [36]. This is significant in our system as the spectral support is wider in the lateral than in the axial direction (in our system respectively 8 µm$^{-1}$ and 0.9 µm$^{-1}$), and since the angular width of the lateral spectrum on the Ewald sphere is large (0.25 rad full angle in our system). The amount of broadening and the resulting shape is, however, also dependent on the lateral frequency content of the object. Figure 7 illustrates this with two extreme examples: a point object, in Fig. 7(a), and a horizontal plane object, in Fig. 7(b). The point object contains many high lateral spatial frequencies, which ISAM maps to lower $k_z$ values. After the lateral inverse DFT, this results in a downshifted and broadened spectrum in the $k_z$ direction. In contrast, a plane object contains a narrow lateral spatial frequency spectrum, see Fig. 7(b). ISAM leaves this axial spectrum unaltered as ISAM remapping occurs at lateral frequencies where there is negligible signal. The dependence of the spectral shape on the local object structure is problematic for SE-OCT as it requires a well-defined spectral shape for spectral normalization. We observed that applying the processing in this order resulted in strong side lobes and a poor axial resolution, especially for structures that had a different spectral content than the used reference spectrum.
Fig. 7. Dependence of the axial object spectrum after ISAM on the object shape. (a) OCT point object intensity before and after ISAM. After ISAM the spectrum bends downward at high lateral frequencies. The spectrum evaluated along $k_z$ is shifted towards lower $k_z$ and broadened. A horizontal plane (b), however, has a narrow lateral spectrum, giving a nearly unaltered axial object spectrum after ISAM. The white dashed lines in the ISAM spectrum indicate the boundaries of the measured OCT data.
The second option does not suffer from axial spectrum ambiguity, but as the signal of small scatterers in the object is smeared out over a large lateral frequency range, the SNR will be low and the spectral estimation not effective. Simulations with this approach showed a lot of artifacts and were thus not satisfactory.
In the proposed method of first applying MIAA followed by ISAM, the data can be uniformly reshaped using the shape of the source spectrum, which is independent of the lateral spectral content of the sample. Even though the axial resolution improvement worsens with the decreasing SNR outside the focal plane, for a large axial range the improvement is still significant, as shown in the current work.
SE-OCT works best in combination with sparse, high-SNR samples [17,18,37]. The performance of the proposed method to non-sparse samples with low SNR will result in less resolution improvement, but will not introduce artefacts or reduce contrast. This was shown with RFIAA in a non-sparse skin sample [18] and retina tissue [38]. Moreover, as MIAA does not add any additional sparsity constraints to RFIAA, MIAA+ISAM is expected to perform similarly in such samples.
A potential disadvantage of MIAA+ISAM is the SNR-caused depth-dependent axial resolution. However, the variation in axial resolution is not as large as the lateral blurring without defocus correction. Moreover, the resolution is still better than the DFT approach over the full depth range.
Using a narrow-band light source greatly simplifies the optical design, especially of the spectrometer. In this work, we used a custom-built OCT setup, with a custom high-resolution spectrometer. However, a spectrometer with one-sixth of the number of pixels (512 instead of 3072 pixels) and a six times lower spectral resolution would be sufficient and not compromise the imaging results, as here only 13% of the maximum imaging depth $z_{\text {max}}$ was used. Moreover, by optimizing the coupling of the fiber source power or increasing the SLD power, the exposure time could be reduced allowing for a faster acquisition with a similar SNR. These adaptions would enhance the applicability of our method.
Another promising hardware adaption could be to add a second SLD and interpolate the signal between the two SLDs using MIAA. If, for example, the blue SLD at 450 nm (EXS210099-03, Exalos) would be added, the 60 nm in between the two spectra could be interpolated and the 70 nm wide spectrum would yield an axial resolution of 1 µm. This would however require a larger spectrometer which complicates the design, like with the green-red combination [14].
The presented results demonstrated that the phase accuracy of IAA-based SE-OCT is sufficient to allow for coherent post-processing using ISAM and CAO. This opens the door to its application to other coherent or phase-based processing methods, such as OCT vibrometry [39], polarization-sensitive OCT, and phase-resolved Doppler OCT imaging at high resolution. CAO could also be combined with induced astigmatism, to improve the SNR outside the focal plane for a better axial resolution without compromising the lateral resolution [9].
The proposed hybrid method of MIAA+ISAM+CAO is especially useful for OCM applications where a high resolution is needed, for example, to reveal sub-cellular details. In the biomedical field, this could be imaging, for example, the cornea [15] or lung tissue that exhibit some sparsity. Outside the medical field, sub-cellular imaging of plants or other small organisms such as fungi could be a useful application area. The strength of MIAA+ISAM+CAO is that the high resolution is obtained with a cost-effective and low-noise SLD source and a simple optic design.
6. Conclusion
In this work, we developed a combined spectral data extrapolation using MIAA with interferometric synthetic aperture microscopy (ISAM) and computational adaptive optics (CAO). We applied this method to data from an OCT setup with a narrow, 6.5 nm, bandwidth SLD source at 510 nm and a high numerical aperture. We obtained a factor 1.8 to 5 axial resolution improvement reaching close to a single µm axial resolution. We obtained a lateral resolution of 0.8 µm over a depth range of $16\times$ (experimental) or $18\times$ (simulations) the nominal depth of field. The close to isotropic resolution resulted in a clear improvement of image quality and optical sectioning on a relevant biological sample. The results show that the accuracy of RFIAA spectral estimation OCT is sufficient for coherent post-processing. Using a computationally efficient implementation of MIAA, the volume image processing time was in the order of minutes. Our findings pave the way for a wider application of cost-effective visible range SLD sources with a narrow bandwidth, as well as for the application of SE-OCT in combination with phase-sensitive OCT imaging.
Funding
Stichting voor de Technische Wetenschappen (16293).
Acknowledgements
This work is dedicated to the memory of Dr. Kostas Angelopoulos who passed away on February 8, 2023. Kostas was a dedicated teaching staff member of the Department of Informatics and Telecommunications and a research member of the Signal and Image Processing Laboratory, of the University of Peloponnese, Greece. We thank Kostas for fruitful scientific discussions and his support in the parallel processing implementation of the spectral estimation algorithms.
Disclosures
The authors declare that there are no conflicts of interest related to this article.
Data availability
Data and code underlying the results presented in this paper are available in Ref. [40].
References
1. Y.-Z. Liu, F. A. South, Y. Xu, P. S. Carney, and S. A. Boppart, “Computational optical coherence tomography,” Biomed. Opt. Express 8(3), 1549–1574 (2017). [CrossRef]
2. J. A. Izatt and M. A. Choma, “Theory of optical coherence tomography,” in Optical coherence tomography, (Springer, 2008), pp. 47–72.
3. A. D. Aguirre and J. G. Fujimoto, “Optical coherence microscopy,” in Optical coherence tomography, (Springer, 2008), pp. 505–542.
4. A. Lichtenegger, J. Gesperger, B. Kiesel, et al., “Revealing brain pathologies with multimodal visible light optical coherence microscopy and fluorescence imaging,” J. Biomed. Opt. 24(06), 1 (2019). [CrossRef]
5. R. Leitgeb, M. Villiger, A. Bachmann, L. Steinmann, and T. Lasser, “Extended focus depth for Fourier domain optical coherence microscopy,” Opt. Lett. 31(16), 2450–2452 (2006). [CrossRef]
6. I. Grulkowski, K. Szulzycki, and M. Wojtkowski, “Microscopic OCT imaging with focus extension by ultrahigh-speed acousto-optic tunable lens and stroboscopic illumination,” Opt. Express 22(26), 31746–31760 (2014). [CrossRef]
7. L. Yu, B. Rao, J. Zhang, J. Su, Q. Wang, S. Guo, and Z. Chen, “Improved lateral resolution in optical coherence tomography by digital focusing using two-dimensional numerical diffraction method,” Opt. Express 15(12), 7634–7641 (2007). [CrossRef]
8. T. S. Ralston, D. L. Marks, P. S. Carney, and S. A. Boppart, “Interferometric synthetic aperture microscopy,” Nat. Phys. 3(2), 129–134 (2007). [CrossRef]
9. S. G. Adie, B. W. Graf, A. Ahmad, P. S. Carney, and S. A. Boppart, “Computational adaptive optics for broadband optical interferometric tomography of biological tissue,” Proc. Natl. Acad. Sci. 109(19), 7175–7180 (2012). [CrossRef]
10. A. Ahmad, N. D. Shemonski, S. G. Adie, H.-S. Kim, W.-M. W. Hwu, P. S. Carney, and S. A. Boppart, “Real-time in vivo computed optical interferometric tomography,” Nat. Photonics 7(6), 444–448 (2013). [CrossRef]
11. A. Lichtenegger, D. J. Harper, M. Augustin, P. Eugui, M. Muck, J. Gesperger, C. K. Hitzenberger, A. Woehrer, and B. Baumann, “Spectroscopic imaging with spectral domain visible light optical coherence microscopy in Alzheimer’s disease brain samples,” Biomed. Opt. Express 8(9), 4007–4025 (2017). [CrossRef]
12. M. Münter, M. Pieper, T. Kohlfaerber, E. Bodenstorfer, M. Ahrens, C. Winter, R. Huber, P. König, G. Hüttmann, and H. Schulz-Hildebrandt, “Microscopic optical coherence tomography (mOCT) at 600 kHz for 4D volumetric imaging and dynamic contrast,” Biomed. Opt. Express 12(10), 6024–6039 (2021). [CrossRef]
13. X. Shu, L. J. Beckmann, and H. F. Zhang, “Visible-light optical coherence tomography: a review,” J. Biomed. Opt. 22(12), 121707 (2017). [CrossRef]
14. A. Lichtenegger, M. Salas, A. Sing, M. Duelk, R. Licandro, J. Gesperger, B. Baumann, W. Drexler, and R. A. Leitgeb, “Reconstruction of visible light optical coherence tomography images retrieved from discontinuous spectral data using a conditional generative adversarial network,” Biomed. Opt. Express 12(11), 6780–6795 (2021). [CrossRef]
15. S. Khan, K. Neuhaus, O. Thaware, S. Ni, M. J. Ju, T. Redd, D. Huang, and Y. Jian, “Corneal imaging with blue-light optical coherence microscopy,” Biomed. Opt. Express 13(9), 5004–5014 (2022). [CrossRef]
16. P. D. Woolliams, R. A. Ferguson, C. Hart, A. Grimwood, and P. H. Tomlins, “Spatially deconvolved optical coherence tomography,” Appl. Opt. 49(11), 2014–2021 (2010). [CrossRef]
17. X. Liu, S. Chen, D. Cui, X. Yu, and L. Liu, “Spectral estimation optical coherence tomography for axial super-resolution,” Opt. Express 23(20), 26521–26532 (2015). [CrossRef]
18. J. De Wit, K. Angelopoulos, J. Kalkman, and G.-O. Glentis, “Fast and accurate spectral-estimation axial super-resolution optical coherence tomography,” Opt. Express 29(24), 39946–39966 (2021). [CrossRef]
19. A. Kumar, W. Drexler, and R. A. Leitgeb, “Subaperture correlation based digital adaptive optics for full field optical coherence tomography,” Opt. Express 21(9), 10850–10866 (2013). [CrossRef]
20. P. Stoica, J. Li, and J. Ling, “Missing data recovery via a nonparametric iterative adaptive approach,” IEEE Signal Processing Letters 16(4), 241–244 (2009). [CrossRef]
21. J. Kalkman, “Fourier-domain optical coherence tomography signal analysis and numerical modeling,” Int. J. Opt. 2017, 1–16 (2017). [CrossRef]
22. T. Yardibi, J. Li, P. Stoica, M. Xue, and A. B. Baggeroer, “Source localization and sensing: A nonparametric iterative adaptive approach based on weighted least squares,” IEEE Trans. Aerosp. Electron. Syst. 46(1), 425–443 (2010). [CrossRef]
23. G.-O. Glentis and A. Jakobsson, “Efficient implementation of iterative adaptive approach spectral estimation techniques,” IEEE Transactions on Signal Processing 59(9), 4154–4167 (2011). [CrossRef]
24. G.-O. Glentis and A. Jakobsson, “Time-recursive iaa spectral estimation,” IEEE Signal Processing Letters 18(2), 111–114 (2010). [CrossRef]
25. T. S. Ralston, D. L. Marks, P. S. Carney, and S. A. Boppart, “Inverse scattering for optical coherence tomography,” J. Opt. Soc. Am. A 23(5), 1027–1037 (2006). [CrossRef]
26. K. C. Zhou, R. Qian, A.-H. Dhalla, S. Farsiu, and J. A. Izatt, “Unified k-space theory of optical coherence tomography,” Adv. Opt. Photonics 13(2), 462–514 (2021). [CrossRef]
27. X. Attendu, R. M. Ruis, C. Boudoux, T. G. Van Leeuwen, and D. J. Faber, “Simple and robust calibration procedure for k-linearization and dispersion compensation in optical coherence tomography,” J. Biomed. Opt. 24(05), 1 (2019). [CrossRef]
28. M. Wojtkowski, V. J. Srinivasan, T. H. Ko, J. G. Fujimoto, A. Kowalczyk, and J. S. Duker, “Ultrahigh-resolution, high-speed, Fourier domain optical coherence tomography and methods for dispersion compensation,” Opt. Express 12(11), 2404–2422 (2004). [CrossRef]
29. N. Nassif, B. Cense, B. Park, M. Pierce, S. Yun, B. Bouma, G. Tearney, T. Chen, and J. De Boer, “In vivo high-resolution video-rate spectral-domain optical coherence tomography of the human retina and optic nerve,” Opt. Express 12(3), 367–376 (2004). [CrossRef]
30. J. de Wit, S. Tonn, G. Van den Ackerveken, and J. Kalkman, “Quantification of plant morphology and leaf thickness with optical coherence tomography,” Appl. Opt. 59(33), 10304–10311 (2020). [CrossRef]
31. G. M. Hale and M. R. Querry, “Optical constants of water in the 200-nm to 200-µm wavelength region,” Appl. Opt. 12(3), 555–563 (1973). [CrossRef]
32. F. Huang, S. L. Schwartz, J. M. Byars, and K. A. Lidke, “Simultaneous multiple-emitter fitting for single molecule super-resolution imaging,” Biomed. Opt. Express 2(5), 1377–1393 (2011). [CrossRef]
33. A. Kumar, T. Kamali, R. Platzer, A. Unterhuber, W. Drexler, and R. A. Leitgeb, “Anisotropic aberration correction using region of interest based digital adaptive optics in Fourier domain OCT,” Biomed. Opt. Express 6(4), 1124–1134 (2015). [CrossRef]
34. L. Glandorf, P.-J. Marchand, T. Lasser, and D. Razansky, “Digital aberration correction enhances field of view in visible-light optical coherence microscopy,” Opt. Lett. 47(19), 5088–5091 (2022). [CrossRef]
35. D. Hillmann, C. Pfäffle, H. Spahr, S. Burhan, L. Kutzner, F. Hilge, and G. Hüttmann, “Computational adaptive optics for optical coherence tomography using multiple randomized subaperture correlations,” Opt. Lett. 44(15), 3905–3908 (2019). [CrossRef]
36. B. J. Davis, S. C. Schlachter, D. L. Marks, T. S. Ralston, S. A. Boppart, and P. S. Carney, “Nonparaxial vector-field modeling of optical coherence tomography and interferometric synthetic aperture microscopy,” J. Opt. Soc. Am. A 24(9), 2527–2542 (2007). [CrossRef]
37. Y. Ling, M. Wang, Y. Gan, X. Yao, L. Schmetterer, C. Zhou, and Y. Su, “Beyond Fourier transform: super-resolving optical coherence tomography,” arXiv, arXiv:2001.03129 (2020). [CrossRef]
38. J. De Wit, K. Angelopoulos, G.-O. Glentis, and J. Kalkman, “Spectral estimation optical coherence tomography for retinal imaging,” in Optical Coherence Tomography, (Optica Publishing Group, 2022), pp. CW1E–5.
39. S. S. Gao, P. D. Raphael, R. Wang, J. Park, A. Xia, B. E. Applegate, and J. S. Oghalai, “In vivo vibrometry inside the apex of the mouse cochlea using spectral domain optical coherence tomography,” Biomed. Opt. Express 4(2), 230–240 (2013). [CrossRef]
40. J. de Wit, G.-O. Glentis, and J. Kalkman, “Computational 3D resolution enhancement for optical coherence tomography with a narrowband visible light source,” Zenodo (2023) https://doi.org/10.5281/zenodo.7870795.
