Computational approach for correcting defocus and suppressing speckle noise in line-field optical coherence tomography images

The trade-off between transverse resolution and depth-of-focus (DOF) typical for optical coherence tomography (OCT) systems based on conventional optics, prevents “single-shot” acquisition of volumetric OCT images with sustained high transverse resolution over the entire imaging depth. Computational approaches for correcting defocus and higher order aberrations in OCT images developed in the past require highly stable phase data, which poses a significant technological challenge. Here, we present an alternative computational approach to sharpening OCT images and reducing speckle noise, based on intensity OCT data. The novel algorithm uses non-local priors to model correlated speckle noise within a maximum a posteriori framework to generate sharp and noise-free images. The performance of the algorithm was tested on images of plant tissue (cucumber) and in-vivo healthy human cornea, acquired with line-field spectral domain OCT (LF-SD-OCT) systems. The novel algorithm effectively suppressed speckle noise and sharpened or recovered morphological features in the OCT images for depths up to 13×DOF (depth-of-focus) relative to the focal plane.


Introduction
Optical coherence tomography (OCT) has been utilized in ophthalmic applications for the past 35 years [1], providing in vivo non-contact volumetric imaging of ocular tissues at high image acquisition rates [2][3][4][5].Monitoring structural changes in the tissue at cellular level caused by various diseases and the cellular response to different therapeutic approaches requires access to high-resolution images.This is particularly important for ophthalmic imaging where morphological changes caused by ocular diseases can eventually lead to blindness.In free space, the OCT axial resolution is determined by the spectral characteristics of the low-coherence light source, while the transverse resolution is determined by the effective numerical aperture (NA) of the imaging lens and the central wavelength and spectral bandwidth of the light source [6].As a result, the OCT axial and transverse resolutions are independent on each other [3].Ultrahighresolution OCT (UHR-OCT) can be realized by employing high-NA optics and broadband light sources [7,8].UHR-OCT allows for in vivo cellular-resolution imaging of ophthalmic tissues, facilitating the study of various eye diseases at their early stages and potentially enabling timely diagnosis [5,9].The use of high NA optics to achieve high transverse OCT resolution results in significant reduction of the depth-of-focus (DOF) due to the inverse square dependence of DOF on the NA.Although this trade-off does not affect the lateral resolution in full-field OCT (FF-OCT) systems, technical challenges associated with FF-OCT prohibit the realization of the axial resolution required for cellular-level imaging.The existing trade-off between transverse resolution and DOF in OCT prevents "single-shot" acquisition of volumetric OCT images with sustained high transverse resolution over the entire imaging depth.
Several solutions have been proposed to tackle the limitations imposed by the trade-off between transverse resolution and DOF in UHR-OCT.In one solution, cellular-resolution OCT data is acquired at multiple depths to cover the region of interest (ROI) in the sample.Subsequently, the volumes are stitched together in the axial direction in post-processing for full-volume results [10].Alternatively, OCT data acquired at different angles with respect to the sample can be integrated in post-processing to synthetize an isotropic point spread function (PSF) that solely relies on the axial resolution of the OCT system, i.e., the coherence length of the light source [11,12].This approach allows for the realization of UHR-OCT with sufficient depth range using only a broadband low-coherence source.However, both solutions require the acquisition of multiple OCT datasets to generate a single full-volume image at cellular resolution.Such approaches are very time consuming, prone to artifact due to incorrect stitching and prevent real-time visualization of the volumetric OCT data.
Another approach is to utilize beam shaping techniques for physical extension of DOF.For example, an axicon lens can be employed to generate a Bessel beam that could be focused to a spot with ∼10 times longer DOF compared to a Gaussian beam [13,14].However, axicon lenses are very lossy and the light intensity along the beam is highly non-uniform which can result in image artifact and loss of signal-to-noise ratio (SNR).Similar DOF extension can be achieved using a phase mask to generate multiple foci in the axial direction [15].However, the limited SNR and the resulting image artifact generated by these methods restricts their applicability for in vivo OCT imaging of the human eye [6].On the other hand, it is possible to counteract defocus as well as higher order aberrations at the time of data acquisition in UHR-OCT using adaptive optics (AO) [12,13].In AO-OCT, total wavefront aberration caused by the optics and the images tissue is measured using a wavefront sensor or a feedback system based on the sharpness of a real-time image.Subsequently, deformable mirrors or other correcting devices such as liquid spatial light modulators or segmented mirrors modify the wavefronts to compensate for the detected aberrations [16].Since aberrations are corrected during data collection in AO, both signal intensity and lateral resolution can be maintained along depth [17].Nonetheless, when the imaging ROI exceeds the DOF of the AO-OCT system, focal plane shifting or computational refocusing in post-processing is required to achieve sharp images of structures at various depths [16,17].Furthermore, AO-OCT systems are bulky, expensive, and difficult to align.As a result, AO-OCT setups are yet to gain widespread adoption in clinical ophthalmic applications due to their expensive and complex design [4].
In contrast to the solutions discussed above, computational methods can be used for the recovery of high-resolution images in UHR-OCT without the requirement for physical modification of the imaging setup [3].In digital refocusing, scalar diffraction models are leveraged for calculating the defocus distances based on a priori information on geometrical properties of the optical design in an OCT system [18,19].Subsequently, the derived defocus terms are employed for designing a 2D phase deconvolution filter in the spatial frequency domain, i.e., Fourier space of enface images.In a similar approach known as interferometric synthetic aperture microscopy (ISAM), OCT data is reconstructed volumetrically by solving an inverse scattering problem in the Fourier space [20].In this method, the 3D spatial frequency data is manipulated to bring all depths into focus, achieving depth-invariant resolution in OCT volumes [3].As a result, ISAM is capable of concentrating the energy distributed in all directions due to defocus and dispersion, and therefore achieves more robust results compared to digital refocusing.Furthermore, digital aberration correction (DAC) can be used for the estimation and correction of total wavefront aberrations in post-processing of complex enface OCT images.Due to their operating principle, DAC methods are regarded as computational equivalents of AO, and therefore, are also referred to as digital adaptive optics (DAO).DAC methods leverage the fact that total aberrations in an enface OCT image can be represented as a multiplicative term in the spatial frequency domain with only a phase component, i.e., unit amplitude.Therefore, given this phase term, aberrations can be canceled through a simple multiplication in the Fourier space.Depending on how the phase aberrations are estimated in DAC methods, they can be implemented in various ways.A group of solutions employ parametrization methods such as Zernike polynomials with variable coefficients to represent the phase aberrations in the Fourier space [4,[21][22][23][24].Next, the coefficients are iteratively optimized using an image sharpness metric such as Shannon entropy.On the other hand, phase aberrations can be non-iteratively estimated using hardware-inspired methods that implement a software version of wavefront sensors employed in AO [6], as described in [25][26][27][28].All the computational methods take advantage of phase information to construct deconvolution filters, and therefore, rely heavily on phase stability of the OCT data.Due to the involuntary eye motion, the applicability of these methods is limited to the data acquired with complex and/or expensive imaging systems that enable phase-stable full-volume OCT measurements in ocular tissues [6,29].Additionally, defocus-induced signal loss restricts the performance of computational methods in recovering high resolution OCT images.This problem has only been mitigated with the aid of AO by maximizing the number of collected photons [17].In addition, the effect of chromatic aberrations in broadband OCT systems is not negligible.Recently, our research group proposed a novel method for correcting chromatic aberrations in broad-band LS-SD-OCT [6] by employing an iterative DAC method for separate aberration correction for monochromatic sub-bands of the broader laser spectrum.However, this method requires using a top surface reflector in contact with the sample for the registration of the sub-band images in post-processing.Such a requirement limits the application of the proposed solution in ophthalmic applications.The drawbacks mentioned above have prevented the wide adoption of current computational methods for recovering high resolution OCT data in clinical ophthalmic applications.As a result, it is increasingly important to address these limitations through a novel computational approach that can be seamlessly integrated into clinical OCT systems for full-volume cellular-resolution imaging of the human eye.
In this paper, we describe a novel computational method for restoring cellular resolution in UHR-OCT images of the human eye and demonstrate its effectiveness by applying it to images of plant tissue (cucumber) and healthy human cornea acquired in vivo with a LF-SD-OCT.

Imaging system
Two generations of an LF-SD-OCT system were used to acquire images from a resolution target, plant tissue (English cucumber) as well as a healthy human cornea, that were used in this study to test the performance of the proposed algorithm.The first generation of the system provided an isotropic resolution of ∼2 µm in free space, ∼18 µm DOF, and a maximum sensitivity of ∼87 dB at ∼2.5 mW optical power and 2,000 fps acquisition rate [30].On the other hand, the second generation LF-SD-OCT system featured an isotropic resolution of ∼2.4 µm, DOF of ∼30 µm, and ∼90 dB maximum sensitivity at ∼2.6 mW optical power and 2,400 fps acquisition rate [31].The volumetric LF-SD-OCT images were comprised of 600 A-scans and 600 B-scans, with ∼1 µm spatial sampling rate, resulting in a ∼600 × 600 µm 2 FOV acquired in ∼0.3 seconds.Imaging human subjects with the LF-SD-OCT system has received full ethics clearance from the University of Waterloo's Research Ethics Board.

Problem formulation
By means of low-coherence interferometry, OCT provides depth-resolved data in various sample types with high axial resolution.However, due to the spatial coherence of the imaging beam and multiple backscattering of light within the samples, the magnitude and phase of the collected signal in OCT will vary randomly.The effect of this randomness in OCT imagery is known as speckle noise.As the size and statistical properties of speckle is a function of the imaging beam spectrum and structural properties of the imaging sample, it can be modeled as a multiplicative noise with an unknown distribution.
Additionally, in order to devise a system-agnostic solution that does not rely on the phase stability of the imaging system, OCT intensity images were used as the inputs to our proposed method.Furthermore, to simplify the complexity of the problem, we chose to work with the enface OCT images corresponding to a single depth layer with an almost uniform PSF across the FOV.As a result of the considerations and simplifications mentioned above, we use the following equation to describe an enface OCT intensity image contaminated with speckle noise: where r is a 2D vector identifying pixels within the images, y(r) is the acquired image, x(r) is the noise-free image, n(r) is the speckle noise, and "•" is the element-wise product operator.
The dependency of the speckle noise on the scattering properties of the samples and the optical design of OCT, and its multiplicative nature, render the problem of recovering x(r) from y(r) especially challenging.Furthermore, image regularization techniques for noise reduction are designed for additive noise models.Therefore, we can convert Eq. ( 1) to an additive equation using a logarithm transformation as following: In the equation above, the "l" subscript denotes representation in logarithm space.Furthermore, the log-scale noise-free image, x l (r), can be represented as the convolution of a blur kernel, h r (r), with the aberration-free image s l (r).This representation models the total nonstationary aberrations in an enface OCT image, hence the "r" subscript to indicate the spatial dependency of h r (r).As a result, the forward model can be formulated as following: where " * " denotes the linear convolution operator.The model in Eq. ( 3) will be used in the following section to estimate the noise-free and sharpened log-scale image s l (r) from the log-scale measurement y l (r).

Correction algorithm
Sharpening defocused OCT intensity images is especially challenging as the speckle noise can result in artifacts that interfere with true structures within the images.Especially, with the goal of obtaining cellular resolution OCT images of the eye for diagnostic purposes, such artifacts can invalidate the results.To address the abovementioned challenges in UHR-OCT images of ophthalmic tissues, we propose a novel approach for maximum a posteriori (MAP) estimation of the noise-free and aberration-reduced image, using a self-similarity-based model as a non-local prior on the speckle noise and the assumption of a Gaussian PSF.A light diffraction model is used for further reconstruction of the noise-free, aberration-reduced image.
In order to account for the spatially-correlated nature of the additive noise term n l (r), we model it as the convolution of a normalized Gaussian kernel, G(r; σ k ), with a 2D i.i.d.Gaussian noise of zero mean and unit variance, g(r) [32], when y l (r) is normalized to the 0 − 1 range.The size of the Gaussian kernel, σ k , depends on the physical properties of the imaging system such as the coherence length of the OCT source.This approach for modeling speckle noise in OCT is similar to that employed in [33], for assuming a stationary model.Moreover, the assumption of stationarity for the additive noise proves to be sufficient [34], when combined with considerations explained later in this section.Furthermore, a uniform (spatially-invariant) PSF was assumed for the sake of simplicity in this work.Consequently, the forward model for the log-scale measurements were modified into the following equation: where the uniform PSF is denoted by h(r).In order obtain a robust estimation of the sharpened and noise-free image, we aim to solve the inverse problem by leveraging non-local priors within a MAP framework as following: ŝl (r), ĥ(r) = arg min where ŝl (r) and ĥ(r) denote the MAP estimations of the sharpened image and the PSF, respectively.Minimizing the first term in Eq. ( 5) pertains to a maximum likelihood (ML) solution without considering the properties of the speckle noise.Simple ML estimation in our case is subject to converging to false local minima leading to unreliable solutions.To address this issue, non-local priors were employed in the form of a regularization term included in the cost function (second term in Eq. ( 5)).Minimizing this regularization term will ensure the structural fidelity of the solution.As such, P n is a prior on the magnitude of the 2D Fourier transformed s n l (r), a patch within the output image.To derive the non-local priors (P n , n = 1, 2, . . ., N p ), we employed a variant of the Block-Matching and 3D-filtering (BM3D) method, tailored for handling correlated noise [32].In this method, the input log-scale image is partitioned into N p patches.Each patch is first 2D Fourier transformed.Subsequently the magnitude of these transform domain images are grouped based on their similarity through a ranking system, described in [32].In the next step, the noise in patches within a subgroup is modeled and mitigated collectively.Leveraging similar patches for deriving noise models allows us to group pixels based on the structural context they belong to, and not necessarily based on their physical proximity.In return, this enables the development of robust priors for removing noise.This approach is akin to the conditional posterior sampling method proposed in [33], where a relevance metric is employed for identifying pixels with similar characteristics.While the design of the OCT systems affects the speckle noise characteristics almost identically across the FOV, patch grouping uniforms the scattering properties of the pixels in the same group.As a result, speckle noise within a subgroup of similar patches can be considered stationary, while the grouping algorithm renders the noise modeling method spatially adaptive.Furthermore, the Gaussian kernel, G(r; σ k ), used in Eq. ( 4) captures the spatial correlation of speckle noise, and is fed into the BM3D algorithm as an input.The size of the Gaussian kernel, characterized through σ k , is mainly determined by the coherence length of the OCT system, reflecting how distant a scattered photon can impact the measured signal at a certain point.Therefore, σ k is expected to be in the ballpark of 1/ (︂ 2 • ∆z where ∆z is the full width at half-maximum (FWHM) of the imaging system's axial PSF.In practice, a near-optimal σ k was found heuristically for a sample enface OCT intensity image based on the improvement in peak SNR (PSNR) and was used for deriving non-local noise priors in all images acquired with the same system.
The processing steps are summarized in Fig. 1.Initially, the raw OCT data goes through standard preprocessing steps for obtaining enface OCT intensity images in log-scale, i.e., y l (r).These steps are shown collectively inside the white dashed box in Fig. 1.Following the preprocessing steps, the non-local priors (P n ) are derived using the method explained above.The speckle model derived using BM3D was then passed to the MAP estimation framework for obtaining a noise-free and sharpened image, ŝl (r).Using the information about the depth location of the input image and optical design of the system, an initial PSF, i.e., h(r), was calculated and passed to the algorithm, which was then iteratively optimized in parallel with the image sharpness.Aside from causing blur, strong defocus in OCT images leads to interfering artifacts due to the introduction of diffraction patterns.In order to address that, the MAP step was followed by diffraction correction using Rayleigh-Sommerfeld (RS) backpropagation method [35], as shown in Fig. 1.To achieve that, the frequency domain representations of the enface images were multiplied by the optical transfer function (OTF) defined in Eq. ( 6).
In the equation above, λ c is the center wavelength of the OCT source, z 0 is the axial distance of the imaged plane from the focal point, F xy is the 2D Fourier transform operator, and (Q x , Q y ) is the pair of spatial frequencies in the transverse directions.Especially for reflecting samples such as the resolution targets, diffraction correction drastically enhanced the details.In practice, the formula in Eq. ( 6) was modified by breaking the symmetry between x and y, as the aberrations in these directions had different characteristics.As such, aberrations in the x direction (width) manifested as uniform defocus while in the y direction (height), diffraction patterns are evident.In order to accommodate for that, the diffraction kernel was stretched along the y direction and then tapered along the x direction using an asymmetrical Gaussian window.

Results
The performance of the novel algorithm was tested on LF-SD-OCT images acquired from a USAF 1951 target, plant tissue (English cucumber) and in-vivo healthy human cornea.Figure 2 shows images of a USAF 1951 resolution target that were acquired using the first generation of the LF-SD-OCT system [30].An in-focus image of the resolution target is shown in Fig. 2(A).The defocused images shown in Fig. 2(B) and Fig. 2(F) were acquired ∼1.8×DOF and ∼5×DOF away from the focal plane, respectively.While the original images were practically speckle-free, we added speckle noise to the images of the resolution target in post processing, using the speckle model described in [34].The image in Fig. 2(F) is more strongly defocused and influenced by diffraction patterns compared to the image in Fig. 2(B), given that it was acquired at an axial distance from the focal plane ∼5 times the DOF of the system (DOF ≈ 18 µm).Moreover, the signal level is evidently lower in this image due to the strong defocus and reduced intensity of the incident optical beam.While signal loss does not pose a limitation in imaging reflective samples such as resolution targets, it can be a significant drawback when imaging scattering  Next, the performance of the algorithm was tested by imaging the soft tissue of an English cucumber near the seed surface.This type of plant tissue has optical transmission and scattering properties similar to corneal tissue and contains cells with sizes ranging from <10 µm to >100 µm.Enface images of English cucumber acquired with the first generation LF-SD-OCT system [30] are shown in Fig. 3. Figure 3(A) shows a photograph of a cucumber slice, while Fig. 3(B) shows a magnified view of inset marked with the red box in Fig. 3(A), corresponding to the soft tissue next to the cucumber seed (red arrow).Compared to the images of the resolution target, speckles are more prominent in the cucumber images, as a result of the presence of scattered light.Additionally, due to the inhomogeneity of the sample, the images are variably defocused across the FOV.Notably, regions with greater defocus exhibit markedly lower SNR due to the reduced incident light intensity.This is best observed in Fig. 3(D) by comparing the bottom-right region of the image with the rest of the FOV.As indicated by the arrows in Fig. 3, the correction algorithm results in cell walls that are well-distinguished from the background and sharpened across various samples.Interestingly, in the areas of extremely low SNR, such as the region in Fig. 3(I) marked with the orange rectangle and orange arrow, the algorithm manages to detect and sharpen cell patterns (corresponding processed images shown in Fig. 3(J)).This is due to the self-similarity-based noise modeling and the fact that in the original images, repeating patterns such the cellular structures provide the BM3D algorithm with a sufficient number of similar patches.Figures 3 (K The algorithm was also tested on healthy human cornea images, acquired in-vivo with the second-generation design of the LF-SD-OCT with a DOF of ∼30 µm in corneal tissues [31].Figure 4(A) shows a cross-sectional image of the human cornea, where the imaging beam was focused roughly in the middle of the epithelium, marked with the red dashed line (depth location Z 0 ).Figures 4(B We have conducted a comparison between our new method and standard image processing algorithms and the results from this study are summarized in Fig. 5.The application of non-local priors within the MAP framework and the RS step have resulted in simultaneous enhancement in the image contrast, sharpness and speckle noise suppression.The "BM3D-MAP" label in Fig. 5 marks the use of BM3D-deriven priors within the MAP framework (rather than sequential application of BM3D and MAP).Full reconstruction (last column in Fig. 5) refers to the use of BM3D-MAP reconstruction method followed by further correction using the RS formula, as demonstrated in Fig. 1.Note that, the resolution target results that lack the RS backpropagation step (Figs.5(B), 5(D), and 5 (G)) contain noticeable diffraction-induced artifacts as it is expected to see only three horizontal bars beside each of the 5 and 6 digits.On the other hand, the artifacts in Figs.5(C) and 5(F) are due to the fact that BM3D denoising was not applied.Additionally, the image in Fig. 5 H is notably sharper than the images in Figs.5(B), 5(C), and 5(E) as a result of using the MAP sharpening step.Finally, using the RS method led to more noticeable improvements in the image quality for the resolution target images compared to the cucumber images (Figs. 4 (G), 4 (H), 4(O), and 4(P)).Similar improvement in the image quality was observed in the cucumber images.

Quantitative analysis
The intensity profiles along the x and y directions in Fig. 3 and Fig. 4 are included to better demonstrate the resolution enhancement and noise suppression abilities of our proposed method.Nonetheless, to provide a quantitative perspective for a more objective evaluation of the performance of the method, we have conducted statistical analysis using the pixel values within the intensity images.As such, Table 1 highlights the SNR and CNR values in nine image pairs from the biological tissues (cucumber and human cornea), obtained using the standard (original) and proposed reconstruction methods.SNR is calculated based on the ratio of average signal level to the standard deviation (STD) of the background noise, whereas CNR captures the ratio  between the contrast (difference between average signal level and the background noise) and the STD of the background noise [36].Notably, the SNR and CNR values are improved in the reconstructed images using our proposed method.a First two samples are selected from the cucumber images.Samples 3 to 9 correspond to the enface images in Fig. 4.

Discussion
The proposed novel algorithm was able to sharpen morphological features in the cucumber tissue such as cell membranes and the grainy structure of the cucumber seed.This led to significantly better visualization of the small cells in the tissue layers closest to the surface of the seed (Fig. 3, regions of interest marked with the orange rectangles and the arrows).
In the images of the human cornea, the algorithm was able to sharpen and improve the contrast of sub-basal corneal nerves (Fig. 4(C)(1)), compared to images acquired with IVCM.This could aid significantly automated morphometry of the corneal nerves using machine learning based image processing algorithms.In the corneal stroma, the novel algorithm was able to correct defocus very effectively in images of stromal keratocytes in a range of depth locations spanning from 1×DOF up to 13×DOF relative to the focal plane (Z0 marked with red dashed line in Fig. 4(A)).Comparison between the original and corrected image at Z7 (Fig. 4(B)(7) and 4C( 7)) demonstrates how effectively the new algorithm can recover morphological details and suppress speckle noise from fairly low SNR images.Comparison of the processed LS-SD-OCT images with IVCM images acquired from similar locations in the healthy human cornea shows excellent correlation in terms of sparseness of the keratocytes distribution within the FOV.Note that size of the keratocytes appears different because the LS-SD-OCT FOV is 600 × 600 µm 2 , while the FOV for the IVCM images is 400 × 400 µm 2 .
It is important to note that since the proposed algorithm utilizes only the OCT intensity data, it is less accurate in recovering out-of-focus-data, especially farther away from the focal plane compared to DAC methods that incorporate phase information [3].This is clear by comparing the images shown in Figs.7), respectively, where some of the keratocytes appear aberrated.In particular, our proposed algorithm was shown effective for restoring images up to ∼13×DOF away from the focal plane, whereas the state-of-the-art DAC methods could achieve diffraction-limited resolution over a depth range of ∼20×DOF [6,37].However, these methods depend on phase-stable OCT measurements in the eye which, in turn, necessitate either a small (FOV) or an advanced OCT setup, such as FF-OCT, that has high phase stability in the XY plane, though not in the XZ direction.Furthermore, the application of the DAC methods in the human eye has been shown only for OCT images of retinal tissues.To our best knowledge, this work provides the first report on recovering defocused OCT images of human cornea acquired in vivo.Nonetheless, the application of our proposed algorithm can be extended to OCT images of retinal tissues as well as other biological samples.
The algorithms in our proposed method were implemented on a desktop computer featuring 64 GB of RAM and an AMD Ryzen 9 7950X CPU.The raw OCT data was preprocessed in MATLAB for generating the enface images.The BM3D algorithm was carried out in Python.Processing each enface image of 600 × 600 pixels took about ∼3 seconds which resulted in a processing time of ∼20 minutes per volume of interest (∼400 slices of 1 µm thickness).The MAP algorithm and the subsequent post-processing steps were carried out in MATLAB.The total run-time of the proposed method is comparable to that of previously proposed DAC methods [4,6].Nonetheless, code optimization, employing a different programming language like C++, or harnessing a machine with more computational power such as a graphics processing unit (GPU) could result in significant improvements in the time complexity of the program.
The processing framework described in this work can be enhanced in terms of output quality in several ways, albeit with an increase in computational complexity.For instance, to enhance the robustness of the algorithm by incorporating volumetric information, the BM3D algorithm can be sequentially applied to enface and B-scan images.On the other hand, while the noise modeling in our proposed method is spatially adaptive and takes into consideration the correlated nature of the speckle, the MAP algorithm was performed assuming a uniform blur across the FOV of the input images, for the sake of simplicity.This assumption is suboptimal because in samples such as corneal tissues, the spatially variant texture and thickness of the sample results in nonstationary aberrations.In order to address this issue, the deconvolution step can be made spatially adaptive by breaking down the input images into patches and sharpening each one independently [38].Additionally, various regularization terms can be used for a MAP estimation of the PSF.Total variation (TV) is an example of a regularization term that emphasizes sharp edges and is suitable for edge preservation and PSF recovery [38][39][40][41].Therefore, TV can be used for obtaining a robust estimate on the PSF and the result can be fed into a MAP deconvolution algorithm with a more suitable regularization term for image recovery such as hyper-Laplacian priors [41].Finally, our method can be used as a complementary tool for conventional DAC methods.While a DAC solution could correct most of the aberrations without a significant modification to the algorithm, our proposed method can address the residual aberrations induced by the broadband source and/or optical misalignments in UHR-OCT systems.

Conclusions
In this work, we described a novel computational approach for restoring cellular resolution in defocused UHR-OCT images of in vivo corneal tissues in the human eye.In the first step, the proposed method leverages BM3D, a nonstationary noise modeling method, for deriving non-local image priors which was then fed into an MAP algorithm for sharpening the images and recovering PSFs iteratively, based on pre-calculated initial PSFs.Subsequently, diffractioninduced artifacts are corrected using a modified Rayleigh-Sommerfeld backpropagation model.The proposed approach was employed to enhance sharpness in enface OCT images in biological tissues including the in vivo human cornea, across a depth range of ∼13×DOF of the imaging system.Given that the proposed algorithm relies solely on intensity information, it is well-suited for OCT applications where inadequate phase stability renders conventional DAC methods ineffective.Therefore, further development of the computational framework introduced in this work can benefit commercial OCT systems for cellular resolution imaging of ocular tissues.Such advancement has the potential to improve diagnostic abilities in clinical ophthalmic applications.

Fig. 1 .
Fig. 1.Block diagram of the full processing steps of the OCT data, including the proposed correction algorithm.The red, yellow, and green boxes indicate the raw OCT data, the preprocessed image input to the correction algorithm, and the final image, respectively.The blue dashed box encompasses the proposed steps for recovering out-of-focus OCT images.
subjects like biological tissues.The corrected images are shown in Fig.2(C) and Fig.2 G, respectively.Figures2(D) and 2(E) show magnified views of the regions of interest (ROIs) in the original and processed images, i.e., Elements 5 and 6 from Group 7, marked with the colored rectangles in Figs.2(B) and 2(C), respectively.Similarly, magnified images in Figs.2 (H) and 2(I) correspond to the ROIs in Figs.2(F) and 2 (G), highlighting Elements 1 and 2 from Group 7 in the original and processed images, respectively.The processed images demonstrate that the proposed algorithm can mitigate diffraction-induced artifacts.

Fig. 2 .
Fig. 2. Resolution enhancement in enface images of a USAF resolution target acquired using a system with DOF ∼18 µm.(A) In-focus image.(B) Image acquired at ∼1.8×DOF away from the focal plane contaminated with speckle noise.(C) The corrected image in (B).(D) and (E) Zoomed-in versions of the regions highlighted by the red boxes in (B) and (C), respectively.(F) Image acquired at ∼5×DOF away from the focal plane contaminated with speckle noise.(G) The corrected image in (F).(H) and (I) Zoomed-in versions of the regions highlighted by the red boxes in (F) and (G), respectively.

Fig. 3 .
Fig. 3. Results of the proposed method in enface images of plant tissues (cucumber).DOF is ∼18 µm.(A) and (B) highlight the area within a cucumber slice that was imaged.(C, E, G, and I) Variably defocused images of cucumber cells.(D, F, H, J) Processed images.The correction algorithm sharpens the cell boundaries even in low-SNR regions.(K) and (L) highlight the ROIs in (G) and (H) marked with red boxes.(M) Intensity profiles along the x direction in (K) and (L).(N) and (O) highlight the ROIs in (G) and (H) marked with red boxes.(P) Intensity profiles along the x direction in (N) and (O).
)(1)(2)(3)(4)(5)(6)(7) show enface LS-SD-OCT images of the corneal stroma at various axial positions with varying levels of defocus and SNR.The corresponding images sharpened with the novel algorithms are shown in Figs.4(C)(1-7), while typical IVCM images from a healthy human cornea are shown in Figs.4(D)(1-7) for comparison.Note, that the LS-SD-OCT FOV is 600 × 600 µm 2 , while the FOV for the IVCM images is 400 × 400 µm 2 .Figure 4(B)(1) shows an enface image acquired at depth Z 1 , roughly 1×DOF (30 µm) away from the focal plane, corresponding approximately to the boundary between the corneal basal cell layer and the Bowman's membrane.Both the OCT and IVCM enface images acquired from this location show presence of sub-basal nerves.The correction algorithm can sharpen the nerves and improve the images contrast (Fig. 4(C)(1)).Enface images obtained from the anterior and middle stroma at depths Z 2 ≈ 2×DOF, Z 3 ≈ 3×DOF, and Z 4 ≈ 5×DOF relative to the focal plane are shown in Figs.4(B)(2-4), respectively.The correction algorithm is able to sharpen images of the stromal keratocytes (bright white spots in the images), suppress effectively speckle noise and improve significantly the image contrast (Figs.4(C)(2-4)).Stromal nerves in Fig. 4(B)(3) are markedly sharpened and brightened in Fig. 4(C)(3) owing to the adaptive noise removal and deconvolution approach.The image shown in Fig. 4(B)(5) was acquired in posterior stroma, ∼7×DOF away from the focal plane.While the structures in this image are sparser than those obtained in anterior stroma, the correction algorithm detects and improves faint and blurred structures (Fig. 4(C)(5)).The enface image shown in Fig. 4(B)(6) corresponds to a depth of ∼9×DOF away from the focus of the imaging beam (∼300 µm deep in the cornea relative to the corneal surface) and shows keratocytes larger in size and more sparsely distributes compared to the anterior stroma.The image shown in Fig. 4(B)(7) was acquired at a depth location ∼13×DOF, corresponding to ∼420 um below the corneal surface.The image shows a very blurred view of the sparsely distributed keratocytes in the posterior stroma.The proposed novel algorithm manages to sharpen the keratocyte images (Fig. 4(C)(7)).Moreover, intensity profiles along the x direction in Figs.4(B)(1-7) and Figs.4(C)(1-7) are shown in Figs.4(X)(1-7).Similarly, Figs.4(Y)(1-7) demonstrate the intensity profile along the y direction in Figs.4(B)(1-7) and Figs.4(C)(1-7).

Fig. 4 .
Fig. 4. Images of the human heathy cornea.(A) A representative B-scan acquired in-vivo with the LS-SD-OCT system shows the layered structure of the cornea.Original and processed enface OCT images corresponding to different depth locations in the cornea are shown in (B1)-(B7) and (C1)-(C7), respectively.For comparison, typical IVCM corneal images of the healthy human cornea are shown in (D1)-(D7).Intensity profiles along the x and y directions in original and processed enface images are shown in (X1)-(X7) and (Y1)-(Y7), respectively.

Fig. 5 .
Fig. 5. Performance comparison between different options for processing enface OCT images and our proposed reconstruction method.(A)-(H) view the Elements 5 and 6 of Group 7 in the USAF 1951 resolution target (top row).Similarly, (I)-(P) refer to a sub-region near the bottom left corner of the sample cucumber image shown in Fig. 3(E) (bottom row).(A) and (I) view the original enface images.(B) and (J) view the results of applying the BM3D method to denoise the input images only.(C) and (K) view the results of only using the RS formula in Eq. (5) to correct for diffraction-induced artifacts.(D) and (L) view the results of using a MAP framework with Gaussian priors to sharpen the input images only.Similarly, images shown in (E) and (M) correspond to the case of using BM3D to denoise the input images and then applying the RS formula.(F) and (N) demonstrate the results of applying the Gaussian-based MAP framework followed by the RS formula.Images in (G) and (O) view the results of using a MAP framework with BM3D-based non-local priors to reconstruct sharp and noise-free images.(H) and (P) correspond to the results of our proposed method 4(C)(6) and 4(C)(7) against those shown in Figs.4(B)(6) and 4(B)(