1. Introduction

Optical coherence tomography (OCT) allows for in-vivo, contactless, non-invasive, volumetric imaging of scattering biological tissues with cellular level resolution and high image acquisition speed. While the axial OCT resolution is defined only by the coherence length of the light source [1], as in other optical imaging techniques, the diffraction-limited OCT lateral resolution is determined by the numerical aperture (NA) of the microscope objective and the beam waist of the optical beam incident on it. The OCT lateral resolution can be enhanced by increasing the effective NA; however, this approach results in a significant reduction in the depth of focus (DOF), which is inversely proportional to the square of the effective NA. Therefore, the optical design of an OCT system has to compromise between the system’s lateral resolution and DOF, depending on the application the system is designed for. In cases that prioritize imaging morphological details of size larger than a single cell, a low NA optical design is used to maintain a relatively uniform lateral resolution over a large imaging depth range [2,3]. In cases where the priority is imaging the cellular structure of biological tissue, high NA imaging optics is typically used [4,5], and depth scanning is required when the regions of interest (ROIs) are larger than the limited DOF. In this case, volumetric OCT images acquired at different depths are merged together using image post-processing algorithms [6].

Different approaches have been developed to overcome the trade-off between lateral resolution and DOF in OCT. The DOF can be physically extended by reshaping the incident plane wave. For example, one method utilizes an axicon lens to generate a Bessel beam illumination that can sustain a uniform lateral OCT resolution over a large depth range [7,8]. Another method uses a phase mask to generate multiple foci along the imaging depth [9]. Both methods can extend the DOF over 10-fold; however, the low SNR prohibits the use of such OCT technology for in-vivo human ocular imaging. An alternative approach to extending the DOF is a Fourier synthesis technique, optical coherence refraction tomography (OCRT) that generates an isotropic synthesized point spread function (PSF) dependent only on the coherence length of the OCT light source, by synthesizing multiple conventional OCT cross-sectional images acquired at different angles [10,11]. However, OCRT requires multi-angle imaging, which limits the image acquisition speed and cannot be directly applied for in-vivo human ocular imaging. Another approach is to digitally refocus OCT images by compensating phase changes in the Fourier space [12–14]. This method requires no modification of the optical design of the conventional OCT technology or the imaging procedure and can achieve spatial invariant lateral resolution over 10-fold DOF in the post-processed OCT images. However, this approach is limited by the lateral phase stability of the OCT system. Therefore, its application for in-vivo OCT imaging, especially in ophthalmology where involuntary eye motion can introduce random phase changes, requires high volumetric acquisition rates, e.g., use of full-field swept-source OCT technology (FF-SS-OCT) [15].

Aside from defocus, higher-order wavefront aberrations caused by imperfections in the optical components and the inhomogeneous structure of the imaged biological tissue, can also deteriorate the lateral OCT resolution. Hardware adaptive optics (HAO) OCT utilizes a deformable mirror to compensate for wavefront aberrations. Aberration-free images of nerve fibers and photoreceptors in the human and animal retina have been acquired in-vivo using this approach [16,17]. However, HAO-OCT systems are associated with significantly higher cost, complex and bulky optical design, and difficulty with the optimal optical alignment. Furthermore, depth scanning or digital refocusing is required for HAO-OCT images when the ROI is larger than its limited DOF. Digital aberration correction (DAC) OCT estimates and corrects the wavefront aberrations in post-processing. Based on how the wavefront aberrations are reconstructed, the DAC algorithms can be categorized into two major groups: iteration-based methods [18–21] and hardware-inspired methods [22–24]. Iteration-based methods estimate the aberration coefficients by optimizing image metric functions, assuming the metric minimum/maximum corresponds to aberration-free images. Hardware-inspired methods, i.e., sub-aperture correlation and digital lateral shearing, are non-iterative and are digital versions of the Shack-Hartmann and lateral shearing wavefront sensors. However, independent speckle patterns lead to inaccurate estimation in the sub-aperture correlation method [25], while the digital lateral shearing method requires single point scatter filtering, which may not be applicable for all types of biological tissues. All DAC methods require high lateral phase stability of the OCT signal, which is dependent on the 3D scanning approach, image acquisition rate, and natural motion of the imaged biological tissue (ocular saccades in the eye, heart rate, breathing rate, etc.) [26,27]. OCT image acquisition rates greater than 2,500 B-scans/second can effectively suppress involuntary eye motion artefacts in ophthalmic OCT images [24]. At such high image acquisition rates, a phase-stable point-scanning OCT device can only image a very small field of view (FOV) or requires an additional motion tracking system [28–30]. Parallel acquisition OCT modalities such as line-scan spectral-domain OCT (LS-SD-OCT) [23] and FF-SS-OCT [15,21], have shown that the use of DAC effectively compensates aberrations in human retinal images acquired in-vivo. However, all published DAC methods designed for OCT are restricted to correction of monochromatic aberrations. Since achieving axial OCT resolution ≤ 2 µm requires light sources with 3 dB spectral bandwidths >150 nm, the chromatic aberrations in such cases are non-negligible. In addition, since the LS-SD-OCT camera sensor records the spatial and spectral information simultaneously, photons with different wavelengths scattered from the same sample position can be projected onto different spatial coordinates of the detector (spatial-spectral crosstalk).

This paper describes a novel approach for digital compensation of spatial-spectral crosstalk and chromatic aberrations in images acquired with a broadband LS-SD-OCT system. The proposed reconstruction method was validated by imaging a standard resolution target, a phantom composed of polystyrene microspheres imbedded in agar gel, and biological tissues such as cucumber and rodent cornea.

2. Factors that degrade the spatial resolution in LS-SD-OCT images

2.1 Virtual interference pinhole effect

In this work, x-direction is defined as the scanning direction of the LS-SD-OCT system, z-direction is the light propagation direction, and y-direction is the direction orthogonal to the xz-plane. As OCT signal is derived from the interference of the electrical fields backscattered from the imaged object (${\vec{E}_s}$) and reflected from the reference mirror (${\vec{E}_r}$), it will reach maximum amplitude only if the two fields are spatially overlapped for all wavelengths of the spectral bandwidth of the optical beam. If the two fields are misaligned in the axial or lateral location (Fig. 1(A) inserts) with respect to each other, the interference fringe’s amplitude drops, which results in a reduced OCT signal. Therefore, the reference electrical field acts as a “virtual pinhole” along the x-direction in LS-SD-OCT. In this study, we named this effect “virtual interference pinhole effect”; however, one should not confuse it with the low spatial coherence illumination [31,32]. Taking the virtual interference pinhole effect into consideration, the effective electrical field backscattered from the sample is (see Supplement 1):

 

Fig. 1. Virtual interference pinhole effect. (A) Schematic of the simplified detection optics of the LS-SD-OCT system. (B) Schematic of the illumination coherent point spread function (cPSF) ${h_i}({x,z} )$ and shifted effective detection cPSF ${h_d}({x,z - Z} )$. (C) Original cPSF ${h_s}({x,z} )$ (left) and refocused cPSF ${\tilde{h}_s}({x,z} )$ (right) when the illumination and effective detection focal planes are overlapped ($Z = 0$). The refocused ${\tilde{h}_s}({x,z} )$ maintains the waist of ${h_s}({x,z} )$ at all depths since phase information is preserved. (D) Original cPSF ${h_s}({x,z} )$ (left) and refocused cPSF ${\tilde{h}_s}({x,z} )$ (right) when the illumination and effective detection focal planes are separated ($Z \ne 0$). ${\tilde{h}_s}({x,z} )$ cannot recover the diffraction-limited resolution in the region between the two focal planes because of phase information loss. F: focal plane, G: grating, L: lens.

Download Full Size | PDF

Here, ${h_{i,d}}$ and $h{\mathrm{^{\prime}}_{i,d}}$ refer to the coherent point spread functions (cPSFs) of the illumination and detection optics in x- and y-direction, $s({x,y,z} )$ is the scattering distribution function of the sample, Z is the separation of the illumination focal plane and conjugate plane of the reference mirror, and ${\otimes} $ denotes convolution.

Since the radius of the illumination beam along the y-direction is much broader than that of the detection beam (420 $\mu $m vs. 1.5 $\mu $m), the effective cPSF of the scattering light can be further simplified as:

It is clear from Eq. (2) that the effective detection focal plane in the x-direction is axially shifted due to the interference effect (Fig. 1(B)), while it is not affected in the y-direction.

The virtual interference pinhole acts as a confocal gate in the x-direction and reduces the necessity of deploying physical apertures; however, it results in phase loss if the focal planes are miss-aligned. Assuming the illumination and detection cPSFs are two identical Gaussian functions with the same beam waist, w, the broadened transverse beam size due to defocusing can be recovered using digital refocusing when the focal planes are perfectly aligned (Fig. 1(C)). However, with a mismatch, the OCT signal’s coherence transfer function at the middle of the illuminating plane and reference conjugate plane is ${H_s}({{k_x},z = Z/2} )\propto wExp({ - k_x^2({{w^2} + {Z^2}/k_0^2{w^2}} )/8} )$, where ${k_x}$ is the spatial frequency and ${k_0}$ is the wavenumber of the central wavelength. For a separation $z = \textrm{Z}/2$, the phase information is completely lost leading to a failure of the digital refocusing approach (Fig. 1(D)).

Alignment of the illumination, detection, and reference focal planes at the same axial position should ideally eliminate phase losses due to misalignment and therefore allow for DAC. However, it is impractical to fulfill this requirement for broadband OCT technology since the axial achromatic aberration is not negligible. Nevertheless, the failure of the DAC to properly compensate monochromatic aberrations is significant only when $|z |> {z_R}$ where ${z_R}$ is the Rayleigh range (See Supplement 1)

2.2 Spatial-spectral crosstalk

In LS-SD-OCT, the spatial distribution of the light backscattered from the imaged object in the y-direction and the spectrum of that light are recorded by an area complementary metal-oxide semiconductor (CMOS) camera. Because of this unique information encoding configuration, it is possible that photons with different wavelengths scattered from the same position $({{y_0}} )$ are projected onto different y-coordinates on the camera: $I({{y_0}} )\to I({y({{y_0},k} ),k} )$, which is defined as spatial-spectral crosstalk in this work.

In LS-SD-OCT technology, the spatial-spectral crosstalk can be caused by the tilted optical axis of the grating with respect to the plane of the camera sensor (Fig. 2(B)), chromatic aberrations in the lateral direction, detection of multiply scattered light, etc. The discussion in this work focuses only on the relative alignment between the diffraction grating and the camera sensor, and the chromatic aberrations. On one hand, the spread of ${y_0}$ to $y({{y_0},k} )$ reduces the lateral resolution in the y-direction (Fig. 2(D)); on the other hand, the distorted spectrum (Fig. 2(C)) results in a broader axial PSF (Fig. 2(E)).

 

Fig. 2. Spatial-spectral crosstalk in LS-SD-OCT. (A) USAF resolution target. (B) Raw data recorded by the camera. The image pattern corresponds to the light reflected by the resolution bars of Group 7 (orange line in A). Tilted lines (yellow box) suggest the existence of spatial-spectral crosstalk. (C) Distorted spectra corresponding to the locations in the camera plane marked with the colored arrows in (B). (D) Enface LS-SD-OCT image of the USAF target. The resolution in the y-direction is affected by the spatial-spectral crosstalk. (E) The axial PSF at the chrome pattern's horizontal edges is much broader than the rest due to distorted spectrum shown in (C).

Download Full Size | PDF

2.3 Chromatic aberrations

Chromatic aberrations in both the axial and lateral directions are not negligible in OCT systems that utilize broadband light sources and cannot be corrected with HAO and monochromatic DAC. Achromatizing lenses have been utilized to correct the axial chromatic aberration in retinal imaging with AO-UHR-OCT [33,34]. However, chromatic aberrations have not been addressed for imaging the anterior human eye (cornea and limbus) with UHR-OCT [3,4,35,36]. Figure 3 shows data acquired by imaging the same region of the USAF resolution target as in Fig. 2, at two different depths separated by ∼15 µm. The spectral lines of the raw camera data corresponding to group 7's horizontal chrome pattern are shown in Fig. 3(A). The axial chromatic focal shift of the two wavelengths, 731 nm and 896 nm, is about 15 $\mu $m which is not negligible since it is twice as large as the system’s Rayleigh range. In addition, the chromatic aberrations in the y-direction lead to a location dependent spatial-spectral crosstalk pattern.

 

Fig. 3. Chromatic aberrations of the LS-SD-OCT in axial (left) and lateral (right) direction. (A) Spectral lines of the resolution target group 7's horizontal chrome patterns. The longer wavelength is focused at z = 0 $\mu $m (A1), and the shorter wavelength is focused at z = 15 $\mu $m (A2). (B) Line profiles of ${I_\lambda }(y )$ of the two wavelengths labeled in (A). (C) Spectral lines of the Group 7 chrome patterns imaged at the side (top) and middle (bottom) FOV in y-direction. The different lateral position's spatial-spectral crosstalk patterns (green and yellow box) are different.

Download Full Size | PDF

The description of the novel digital approach to correcting for spatial-spectral crosstalk and chromatic aberrations in LS-SD-OCT images, as well as the images of phantoms and different biological tissues are presented in the sections below.

3. Methods

3.1 Broadband LS-SD-OCT system

The optical design of the LS-SD-OCT system used to collect the data presented in this manuscript was similar to that of a recent publication [37]. Briefly, the LS-SD-OCT system is powered by a supercontinuum laser (SuperK, NKT Photonics) and data is collected with an area camera (Dimax S4, PCO). The system provides 1.7 µm axial resolutions in free space and sensitivity of ∼ 89 dB near the zero-delay line with ∼10 dB SNR roll-off over scanning range of 1 mm for ∼2.5 mW imaging power. The system’s phase stability is about 0.01 rad (std) at 36 dB. The system’s design and scanning protocol allows for acquisition of a single 3D tomogram with a lateral field-of-view (FOV) of 0.86 mm ${\times} $ 0.56 mm (x ${\times} $ y) and 1.27 mm (z) depth scanning range in free space within 0.17s. All 3D images included in this manuscript consist of 500 B-scans.

3.2 Sample preparation

The proposed method for correction of the spatial-spectral crosstalk and chromatic aberrations was tested on images of an USAF resolution target, a microsphere phantom, as well as biological tissues such as cucumber and a post-mortem rat cornea, acquired with the LS-SD-OCT. The phantom was prepared by imbedding polystyrene microspheres into agar-based gel. The agar gel was prepared by mixing water and agar powder (Sigma Aldrich, St. Louis, US) at a weight ratio of 20:1. The agar powder was first dissolved in water at temperature 80 ∼ 90 ${^ \circ }C$ and the agar gel was cooled down to 70 ${^ \circ }C$ before dispersing polystyrene microspheres with a diameter of 3 $\mu $m (Polybead, Polysciences, Warrington, US).

A Spague-Dawley female rat was euthanized, and volumetric images of its cornea were acquired with the LS-SD-OCT system immediately after euthanasia. Ex-vivo imaging of the rat cornea was necessary due to COVID-19 related restrictions that prohibited research studies on live animals at the time these experiments were carried out. The post-mortem animal imaging procedure was approved by the University of Waterloo Office of Research Ethics. Glass coverslips were placed on the surface of the microspheres’ phantom, cucumber, and the rat cornea during the imaging procedure in order to reduce specular reflections from the surface of the imaged object and to allow for axial phase correction and registration.

3.3 Iteration-based DAC approach

The iteration-based DAC is an optimization problem with the goal of finding the aberration coefficients at which a chosen image quality matric reaches its global minimum/maximum. The OCT enface signal corrected by a guessed wavefront error is:

An appropriate image quality metric is essential to the optimization process, and its global minimum/maximum should correspond to the fully aberration-corrected image. In this study, the image metric functions of the scattering [38] and reflecting OCT signals are chosen as:

The intensity based matric S_scat is commonly used in the DAC methods for scattering biological tissues. However, it fails to represent the image quality reflecting signals (e.g., resolution target) as its global minimum is far away from true aberration coefficients. In contrast, the probability based Shannon entropy S_ref represents correctly the image quality of the reflecting signal but is much noisier than S_scat in the scattering case (see Supplement 1). We limit the use of S_ref/S_scat in correcting the monochromatic aberrations of the resolution target/microsphere phantom and biological tissues.

In this study we choose cuckoo search as the iteration algorithm [39]. Cuckoo search is a typical swarm intelligent algorithm which uses multiple agents to collectively explore the parameter space. The update rule of cuckoo search has two phases: a global random walk and a local random walk. The global random walk combines stochastic Lévy flights and a deterministic attraction to the global best location $\vec{g}\ast $:

As a stochastic algorithm, cuckoo search is more robust to find the global optima of the image metric than deterministic algorithms (e.g., gradient-based algorithms), especially in low SNR regime where there are numerous local minima/maxima (see Supplement 1). Therefore, we choose cuckoo search to make sure the determined aberration coefficients reliable and robust, at the cost of relatively slow converging rate.

3.4 OCT signal reconstruction procedure

The virtual interference pinhole effect described in section 2.1 requires better optical alignment of the LS-SD-OCT system. It cannot be corrected digitally because of the phase information loss. However, spatial-spectral crosstalk and chromatic aberrations can be digitally compensated by independently correcting the monochromatic aberrations in the sub-band and sub-region tomograms and registering them afterward. Here, “sub-band” refers to narrower spectra, and “sub-region” refers to a reduced FOV in the xy-plane.

The flowchart in Fig. 4 illustrates detailed the OCT signal reconstruction procedure in this study. The interference fringes of three sub-bands, ${I_{A,B,C}}({x,y,k} )$, were first obtained by applying three Gaussian filters on the background-subtracted and spectrum resampled signal. Each sub-band fringe was dispersion compensated, and then Fourier transformed to generate the OCT signal $I{\mathrm{^{\prime}}_{A,B,C}}({x,y,z} )$. The weighted average of phase difference between adjacent B-scans due to sub-micrometer involuntary axial motion was calculated as [23,29]:

For each sub-band and sub-region OCT image, the aberration coefficients of the wavefront error at every 5 depth locations were estimated with an iteration-based DAC method. At each depth, the enface image was first coherently averaged with 4 adjacent enface images:

The mean enface signal was then transformed into the Fourier space. Next, the spatial frequency representation was multiplied by an aberration phase mask, and inversely Fourier transformed back to its spatial representation. The iteration was repeated until the image metric stopped improving or the iteration number reached its maximum value. The estimated aberration coefficients were then smoothed along depth using the robust locally estimated scatterplot smoothing (RLOESS) method [40]. Next, the aberration-free tomograms were generated as:

For each sub-region, the depth-dependent lateral shifts of sub-band A and C with respect to sub-band B were determined for every 5^th depth pixel using a subpixel image registration algorithm [41,42]. Once the sub-pixel lateral shifts $\tilde{s}_{x,y}^{A,C,m}(z )$ were smoothed along the axial direction, the sub-band A and C were digitally shifted and registered with sub-band B:

Once the lateral registration was done, the corrected interference fringes were generated via inverse Fourier transform of the three sub-band tomograms along the z-direction:

Because of the inevitable residuals of the dispersion, the phase of the corrected inference fringes of the three sub-bands are discontinuous and proper phase compensation is required [17,43]. To determine the phase discontinuities, the interference fringe of the top air-glass interface was isolated:

Finally, the corrected OCT tomogram was generated by Fourier transforming the phase registered full interference fringes.

4. Results

The results generated with different OCT image reconstruction methods are compared in this section. We have assigned the following names to the different methods:

 

Fig. 4. Flowchart of the OCT signal reconstruction procedure.

Download Full Size | PDF

4.1 3D resolution recovery in resolution target and microspheres phantom images

The recovery of the diffraction-limited lateral resolutions and ultra-high axial resolution in the LS-SD-OCT images was first demonstrated in an USAF resolution target and a microspheres phantom. Metric ${S_{\textrm{ref}}}$ was used in the DAC procedure when correcting the aberrations of the resolution target image. Only the sub-region containing groups 6 and 7 patterns is presented in Fig. 5. The enface and B-scan images of the original, full-spectrum corrected, sub-bands corrected and lateral registered without and with axial phase registration OCT tomograms of the resolution target are shown in Fig. 5(A-D). The aberration coefficients of the enface image with the maximum intensity were determined to the 5th order and were then applied to the whole 3D tomogram. The enface image shown in the figure corresponds to the imaging depth at which the OCT signal reached its peak, and the B-scan images correspond to Group 7's horizontal patterns indicated by the black line ( Fig. 5(A)).

 

Fig. 5. 3D recovery of the OCT diffraction-limited lateral resolution and coherence length confined axial resolution in the resolution target image acquired with the LS-SD-OCT system. (A-D) Log-scale enface and B-scan images of the original, full-spectrum corrected, sub-bands corrected and registered without and with axial phase registration OCT tomograms of the resolution target. The B-scan was acquired at the location of Group 7's horizontal patterns indicated by the black line in (A). (E) The real and imaginary part of the interference fringes at Group 7 Element 6, indicated by the yellow arrow in (D). (F) Phase registered interference fringes at the same location of (E). (G) Axial PSF at the location labeled in (D).

Download Full Size | PDF

The spatial-spectral crosstalk impacts the full-spectrum corrected result in both the lateral and axial directions. In the y-direction, the impact is shown as the blurred horizontal patterns in the enface image. In z-direction, the impact is shown as the broadened signal distribution and the lateral shifts of the signal pattern at different depths in the B-scan image. The broadening of the signal corresponds to the distorted spectrum, and the lateral shifts might be generated by changes in the optical path length when the light traveled to neighboring pixels in the y-direction. Direct summation of the three laterally registered sub-band tomograms results in a reduced depth resolution due to the narrower bandwidth of each band and axial shifts of the A-scan if the dispersion is not fully compensated. With proper phase registration (Fig. 5(E) and 5(F)), the ultra-high axial resolution can be recovered. The intensity of the A-scan signals, $\hat{I}(z )= {|{{E_r}{E_s}(z )} |^2}$, from the chroma pattern of group 7 element 6 is compared in Fig. 5(G). The proposed method achieved similar FWMH with the original signal, with a much higher peak intensity.

Metric ${S_{\textrm{scat}}}$ was used in the DAC procedure when correcting the aberrations of the microspheres image. Typical B-scans of the original, full-spectrum corrected, and fully corrected OCT images are compared in Fig. 6(A), and a selected ROI near the focal plane was zoomed for better visualization in Fig. 6(B). The axial PSFs of the microsphere labeled by the green line are shown in Fig. 6(C). The tilted cross-section of the microspheres away from the focal plane (Fig. 6(B1)) shows the spatial-spectral crosstalk artifact again on the OCT image quality. The full-spectrum correction improved the lateral resolution to some extent, while its axial resolution is much larger than the coherence gate of the light. Our proposed method achieved the best image quality with the highest spatial resolutions. The proposed method produced noticeable sidelobes in the axial PSF, which resulted in ghost pairs around the true microspheres (Fig. 6(B3)).

 

Fig. 6. 3D resolution recovery in the microspheres phantom images. (A) Log-scale B-scan image of the original, full-spectrum corrected, and sub-bands corrected and fully registered OCT images. (B) Enlarged microspheres images of the ROIs labeled in (A). (C) The intensity of the axial PSF scattered by the microsphere labeled by the green line in (B1). (D) Enface image of the original, full-spectrum corrected, and sub-bands corrected and fully registered OCT images at the depth (z = 360 pixel) labeled by the yellow dash line in (A1). (E) Zoomed ROIs selected in (D). (F) Amplitude of OCT signals in x- and y-directions scattered by the microsphere labeled (green solid and dash lines) in (E1). (G) Lateral shifts between sub-band A and B at different sub-regions labeled (red and blue boxes) in (D1). The labeled number indicates the index of the sub-region. (H) The refocusing coefficients, ${a_{2,2}}$ and ${a_{2,0}}$, of the three sub-bands at sub-region 10. Both the enface and B-scan images were incoherently averaged with 2 adjacent neighborhoods. The lateral shifts and refocusing coefficients were smoothed along depth. Scale bars are 100 $\mu $m (original size) and 50 $\mu $m (zoomed).

Download Full Size | PDF

The enface images of the three tomograms at the labeled depth (Fig. 6(A1), yellow dash line) are shown in Fig. 6(D), and the selected ROIs were enlarged and shown in Fig. 6(E). The full spectrum corrected image showed a decent image quality improvement (Fig. 6(D2)), however, still not as good as the proposed method (Fig. 6(D3)). In addition, due to the spatial-spectral crosstalk, a single microsphere can be split into two in the full-spectrum corrected images, as indicated by the yellow arrow, while the proposed method showed the true distribution. The lateral cPSFs of the labeled microsphere (Fig. 6(E1), green lines) along the x- and y-directions were shown in Fig. 6(F). The cPSF generated by the proposed method shows the narrowest waist in both lateral directions. The depth-dependent lateral shifts between sub-band A and B ($\tilde{s}_{x,y}^{A,z,m}$) at different sub-regions are shown in Fig. 6(G). The y-shifts are from the 2nd column sub-region 2, 6, 10,14 (Fig. 6(D1), red boxes); while the x-shifts are from the 3rd row sub-region 9, 10, 11,12 (Fig. 6(D1), blue boxes). The y-shifts of the 4 sub-regions are well separated from each other across the whole depth range, and each y-shift's variance is small compared with the gaps. The x-shifts are about 0 at the depth far away from the focal plane; however, there are significant changes near the focal plane. The lateral shifts between images of sub-band C and B ($\tilde{s}_{x,y}^{C,z,m}$) are similar to $\tilde{s}_{x,y}^{A,z,m}$. The x- and y-shift curves suggest that the sub-region division and frame-by-frame lateral shift registration are essential for reconstruction of the OCT signal in the broadband LS-SD-OCT system. Figure 6(H) shows the refocusing coefficients, ${a_{2,2}}$ and ${a_{2,0}}$, of the three sub-bands at sub-region 10. The obvious differences between the refocusing coefficients of three sub-bands reinforce the importance of applying DAC within narrow sub-bands instead of the full spectrum. As expected, ${a_{2,0}}$ (refocusing coefficient in the y-direction) are linear to the depth z. Because the system's cPSF along y-direction only depends on the detection optics, and its defocusing phasor is proportional to $z\lambda $. In contrast, the calculated ${a_{2,2}}$ (refocusing coefficient in the x-direction) is nonlinear near the focus plane, especially in sub-band C, because of the illumination/effective detection plane mismatch.

4.2 Ex-vivo imaging of biological tissues

The proposed correction approach was also tested on biological tissues (cucumber and rat cornea). Metric ${S_{\textrm{scat}}}$ was used in the DAC procedure. To balance the processing time and image quality, the aberrations were corrected to the 2nd order in DAC. The LS-SD-OCT system is able to resolve the small immature cells located at the surface of the cucumber seed. The diameter of these immature cells is typically about 10 $\mu $m, which is smaller than the epithelium cells of the human cornea. However, due to the limited DOF and the orientation of curved surface of the cucumber seeds, only a small fraction of the immature cells can be resolved (Fig. 7(A), top). The proposed reconstruction method has successfully refocused the images, and the small cells in the regions indicated by the red arrows and yellow box can be resolved (Fig. 7(A), bottom). Figure 7(B) shows a slightly tilted enface projection of another cucumber seed with a different orientation. Although the large cells can be resolved in the original image, the corrected one shows thinner cellular membranes (yellow arrows). In addition, the small immature cells are better visualized in the corrected image, as shown in the zoomed ROI.

 

Fig. 7. Immature cucumber cells of the original (top) and corrected (bottom) OCT images. (A) Tilted projection images with the zoomed ROIs (yellow box). (B) Tilted projection images of a different seed and the enlarged ROIs (red box). All images were generated with Amira software. Each image was incoherently averaged with two adjacent neighborhoods. Scale bars are 100 $\mu $m (original size) and 50 $\mu $m (zoomed).

Download Full Size | PDF

In the case of ex-vivo rat cornea imaging, the imaging focal plane was first positioned at a depth near the boundary between the epithelium layer and the Bowman's membrane. All the five layers of the rat cornea are resolved in the original B-scans (Fig. 8(A), top). The image quality improvement in the corrected B-scan (Fig. 8(A), bottom) is not as noticeable as in the microspheres’ phantom images and the cucumber seeds images because of the lack of small morphological structures. Nevertheless, the reduced cross-section of the stromal nerve (yellow arrows) located in the middle of the stroma is clear evidence of the effectiveness of the proposed reconstruction method. In the enface images, the superiority of the corrected image is more apparent. The local contrast has been enchanted in the enface image near the focal plane (Fig. 8(B)), while the resolution is recovered in the out-of-focus enface image (Fig. 8(C)).

 

Fig. 8. Original (top) and corrected (bottom) OCT images of ex-vivo rat cornea. (A) A typical B-scan image of the rat cornea. (B) The enface image underneath the boundary of the epithelium layer and Bowman's membrane. (C) Enface image in the middle stroma. All images were generated with Amira software. Each image was incoherently averaged with two adjacent neighborhoods. Ep: epithelium; BM: Bowman's membrane; St: stroma; DM: Descemet's membrane; Ed: endothelium layer. Scale bars are 100 $\mu $m.

Download Full Size | PDF

Figure 9 shows maximum intensity projection images of the rat corneal endothelium layer. The focus position was shifted to the posterior end of the cornea. However, due to the tilt and curvature of the endothelium layer, only a small fraction of the cells are in focus. The images of the endothelium layer was first corrected based on the extrapolation of the aberration coefficients determined within the stroma, followed by a flattening algorithm that shifted the endothelium surface to the same depth with the sub-pixel resolution. Both the high and low frequencies in the 3${\times} $3 sub-regions were then filtered out with a bandpass filter, and the contrast of the 9 sub-regions was adjusted to the same level. The reconstructed images reveal the boundaries of the endothelial cells, which are blurred or missing in the original images.

 

Fig. 9. Flattened, enface LS-SD-OCT images of a rat’s corneal endothelium. Original image (A) and corrected (right) enface maximum intensity projection (MIP) of the endothelium layer. (C-F) Enlarged endothelium images of the selected ROIs labelled in (A) and (B). Scale bars are 100 $\mu $m (original size) and 50 $\mu $m (zoomed).

Download Full Size | PDF

5. Discussion

This study investigated some factors that cause degradation in the spatial resolution of volumetric images acquired with a broadband LS-SD-OCT system and impede the success of DAC algorithms to restoring the diffraction-limited lateral resolution. These factors are associated with the system's optical design, alignment, and the broadband nature of the low coherent light source necessary to achieve high axial OCT resolution. The virtual interference pinhole model predicts the phase information loss in the x-direction, which is unique among different configurations of OCT. In point-scan OCT systems, the fiber system separates the light collection process from the interference process. Its cPSF is determined by the optics in the sample arm, and the symmetry of the illumination and detection arms guarantees no mismatch between the two focal planes. In a free-space LS-SS-OCT design, the spectral information is encoded in time instead of space. Therefore, it is the pixel size of the linear camera that confines the region at which the interference signal is collected. FF-SS-OCT deploys a wide-field configuration in both lateral directions; thus, its effective cPSF is the detecting cPSF itself. Some FF-SS-OCT use low spatial coherence light source to reduce the spatial crosstalk [31,32], and the reduced coherence area sometimes is also called “virtual pinhole”. However, the two “virtual pinhole” effects are different. In the broadband LS-SD-OCT system, it is impossible to align the illumination, detection, and reference focal planes at the same axial position for all wavelengths, considering the axial chromatic aberrations of the three paths can be different from each other. In practice, if the separation is within the Rayleigh range, the resolution loss in x-direction due to phase loss is not evident.

The spatial-spectral crosstalk results directly from the mismatch between the optical axis orientations of the camera detection array and the transmission grating. The angle between the two axes in our LS-SD-OCT system is typically less than 7.2 arc minutes but is not negligible. Due to the mechanical drift of the camera, perfect alignment is unrealistic. Moreover, the crosstalk curves are location-dependent because of the chromatic aberration in y-direction; thus, pre-calibrating the crosstalk before each imaging session is strenuous and impractical.

The chromatic aberrations have not been widely addressed for OCT imaging of the anterior segment of the eye or other biological tissues. The distortion caused by chromatic aberrations might be insignificant in OCT systems powered by narrowband light sources, especially in morphological studies where lateral resolution greater than 10 $\mu $m. However, the chromatic aberrations cannot be ignored in the case of ultrahigh resolution LS-SD-OCT system which offers nearly isotropic ∼2 µm resolution in all three dimensions. In our system, the axial chromatic aberration induced focal shift between the maximum and minimum wavelengths is larger than 15 $\mu $m, about twice the Rayleigh length. The chromatic aberrations in x- and y-directions behave differently; the former is depth dependent while the latter is not. The maximum lateral shifts induced by chromatic aberration were about 3 $\mu $m and 2 $\mu $m in x- and y-direction respectively, comparable with the diffraction-limited resolutions. A custom achromatizing lens has been used in retina imaging with point-scan AO-OCT to compensate for the axial chromatic aberration [33,34]. A similar design can be applied in our system; however, unlike the dual-pass configuration in the point-scan system, the illumination and detection optics are not symmetric in the LS system, and their axial chromatic aberrations may need to be compensated separately by different achromatizing triplets. By far, there are no simple optics elements to correct for the chromatic aberrations in lateral directions; thus, digital alignment is necessary to compensate for the lateral shifts [44].

Spatial-spectral crosstalk and chromatic aberrations can be reduced by independently registering the aberration-corrected sub-band OCT images in different sub-regions. Anisotropic aberrations across the lateral FOV have been observed and separately corrected in different studies [21,45,46]; however, these studies have been limited in the monochromatic regime. We take a further step to split the full spectrum into different sub-bands. The lower limit of the sub-bands’ number depends on the special-spectral crosstalk. The proposed method assumes that there is no obvious crosstalk within each sub-band. In our system the spatial-spectral crosstalk spanned across 4-5 pixels in the y-direction. We chose 3 bands in this study. While increasing the sub-bands’ number may further reduce the crosstalk and chromatic aberrations, the benefit is not significant considering the requirement for additional computational time. In addition, there is an upper limit on how many spectra can be used because the accuracy of the DAC method is related to the SNR of the sub-band images. Although, the SNR of the image can be increased by averaging multiple volumes, the lateral phase stability will most likely be not preserved during in-vivo acquisition of multiple volumes.

Since the defocusing phase is proportional to the wavelength, a fully defocus correction should be wavelength dependent. To our best knowledge, Interferometric Synthetic Aperture Microscopy (ISAM) is the only method by far to account for that Ewald sphere curvature [13]. Nevertheless, a comparison study showed no significant difference between the improvement generated by ISAM and monochromatic refocusing algorithms in a narrowband OCT system [47]. Thus, we treat the defocus aberration as monochromatic in each sub-band. Theoretical analysis shows that the LS OCT's 3D coherence transfer function (CTF) is not separable in the two transverse spatial frequencies [48]. However, since our system's effective NA is relatively small (∼0.13), we assume this limitation does not critically affect our analysis and algorithm in each sub-band.

One problem that prevents generating an aberration-free image via metric optimization is the unpredictable local minima. Gradient-based algorithms are intrinsically easy to get trapped in local minima. The SI-based algorithms, particularly the cuckoo search method, show a better convergence to the global minimum at the cost of increased function evaluation numbers. It took about 3 hours to generate the final corrected OCT volume using a customized MATLAB program with a MacBook Air (M1, 16 GB memory). No parallel programming or 3^rd party library were used. The time can be reduced by using fewer enface layers for determining the aberration coefficients and more powerful computers.

Incompletely compensated dispersion results in the axial shifts of the three sub-band tomograms. Phase registration is necessary to correct the phase discontinuity between sub-bands. We estimated the phase differences based on the interference fringe of the air-glass interface, assuming the dispersion is depth invariant and the phase discontinuity is constant across the whole depth range. For future in-vivo applications where coverslip cannot be applied, the phase differences can be estimated based on the axial shifts between the sub-band’s tomograms [17,43].

Like most DAC methods, maintaining lateral phase stability is crucial for the success of our proposed image reconstruction method. In ex-vivo validation, the motion is within sub-micrometer level, allowing for straightforward correction of phase shift between adjacent B-scans. However, in the case of in-vivo ophthalmic imaging, involuntary eye movement may result in larger phase errors and potential inaccuracies in the phase correction. Our LS-SD-OCT system requires approximately 0.2 seconds to acquire a single volume, and it is not guaranteed that motion will be completely absent within such timeframe, especially when averaging of multiple volumes is required for resolving low contrast morphological details. Effective motion and phase error correction methods still need to be developed before implementing our method for in-vivo applications.

6. Conclusions

In summary, we demonstrated an OCT reconstruction method that can restore the diffraction/coherence limited resolutions in a broadband LS-SD-OCT system. An iteration-based DAC algorithm was independently applied to the sub-band tomograms at different sub-regions to correct monochromatic aberrations. The spatial-spectral crosstalk and chromatics aberrations were effectively suppressed by registering the sub-band tomograms at sub-regions. In the microspheres phantom, the proposed reconstruction method showed superior resolutions compared with the full-spectrum monochromatic DAC method. Diffraction-limited transverse resolutions were achieved over a depth range over 300 $\mu $m, which is about 20${\times} $ larger than the DOF. 3D reconstruction of biological tissue images demonstrated the ability to recover the scattering contrast of cells and tissues, which were indistinguishable in the original tomograms. The reconstruction method for the broadband LS-SD-OCT enables fast 3D cell and tissue with isotropic micrometer level resolution across a large depth range and makes it a suitable tool for functional imaging, flow measurement, and dynamic signals with proper motion and phase error correction methods.