Optical coherence tomography (OCT) is a powerful technique that provides non-invasive images with micrometer resolution and millimeter depth. To assess the tissue microstructures at the cellular and subcellular level, advancements^1-6 have been made to improve spatial resolutions and image contrast. However, high-resolution OCT images often suffer from a broaden or distorted point spread function (PSF) due to the inherent trade-off between resolution and depth-of-focus, artifacts resulting from strong scattering, optical component deficiencies, speckle noise from coherent light source, and imaging system's noise. Therefore, enhancing image quality of the OCT system holds significant importance. Without introducing more optical components, many software-based approaches have been harnessed to recover the degraded OCT images by deterministic models^{7, 8} or a class of evaluation metrics^9-11.

Deconvolution, an algorithm extensively used in fluorescence microscopy, has shown promising potential in mitigating transverse and axial resolution blurring in OCT images. Numerous attempts have been made to apply it in OCT, leading to enhanced image quality in certain scenarios^12-18. However, when amplitude-based data is used to invert a coherent process, the image quality will be inevitably deteriorated as the operation creates sidelobes from the image noise around the scatterers. Moreover, closely adjacent scatterers in en-face images are difficult to resolve, limiting the use of deconvolution for OCT imaging scattering tissues¹⁹. Transverse and axial resolutions are separable in OCT systems, where the axial resolution is governed by the central wavelength and the bandwidth of the detected light while the transverse resolution is limited by the effective numerical aperture. With a separable one-dimensional (1-D) scheme, each depth slice or A-line of an OCT frame is deconvolved by the unregularized Richardson-Lucy algorithm. Some regularization constraints, such as Tikhonov-Miller regularization²⁰ and Total Variation regularization²¹, have been used to suppress noise amplification, yet these methods require prior knowledge of the measured data and can result in either smooth edges or a contrast loss in the image restoration.

In this work, we propose a post-processing algorithm, Deconv-RPM. The Fourier or pupil plane of a deconvolved frame is separately multiplied with multiple random phase masks before the products are inverse Fourier transformed to the image plane and summed together. Random phase masks (RPM) are employed to remove the artifacts while preserving the resolved features with deconvolution, a concept vaguely reminiscent of numerical random masks methods applied in digital holography for image denoising^22-24. Our approach, verified through simulation testing, has achieved a 2.5-fold reduction in the full width at half-maximum (FWHM). When applied to a turbid sample, Deconv-RPM exhibited unprecedented performance in cross-sectional imaging compared to the unprocessed images. Deconv-RPM was used to isolate scattering signals from nonkeratinized epithelial cells exvivo and blood cells invivo, using both standard and high-speed line scan cameras. This robust and reliable method is expected to pioneer new possibilities for OCT in biomedical applications.

In this section, we elaborate on the Deconv-RPM method that is schematically illustrated in Fig. 1. For coherent imaging, image formation of one B-scan in OCT can be described as a convolution of the electric field of the original undistorted sample object u_obj(x, z) and the 2-D PSF h(x, z) followed by adding noise n,

1uimg(x,z)=uobj(x,z)⊗h(x,z)+n,

where ⊗ denotes the convolution operation.

The backscattered magnitude of the conventional OCT frame |u_img(x, z)| is obtained by the Fourier transform of the raw interferogram in k_z-space. As demonstrated, the convolution kernels can be approximated by absolute values instead of the complex values^{14, 25}. The OCT PSFs are separable in transversal and longitudinal coordinates. That is h(x, y, z) = h(x)·h(y)·h(z), and can be approximated by

2h(x,y,z)=exp[−x2wr(z)2]⋅exp[−y2wr(z)2]⋅exp[−4ln2(z/lC)2],

where wr(z) is the spot radius of the focused beam and lC is the coherence length. As illustrated in the flowchart, the Deconv-RPM method primarily consists of two processing steps: the first step involves deblurring through deconvolution, and the second step improves image quality using numerical random phase masks.

In step 1, Richardson-Lucy deconvolution algorithm is implemented on one B-scan |u_img(r)|, performed by the built-in MATLAB function [deconvlucy ()],

3ot+1(r)=ot(r)⋅(|uimg(r)|ot(r)⊗|h^(r)|⊗|h^T(r)|),

where the pixel index r represents any one of the coordinates x or z for B-scans and x, y or z for C-scan, ot(r) is the estimated object image after t-th iteration. |h^(r)| is the theoretical 1-D absolute value of the PSF, and |h^T(y)| is its transposed matrix. In the Richardson-Lucy deconvolution method, |h^(r)| is typically set as space-invariant. The primary challenges for deconvolution methods to resolve images in OCT lie in the absence of an accurate complex PSF due to experimental constraints and the inevitable ringing artifacts that occur during the deconvolution process.

In step 2, the deconvolved image o(r) is further modulated by the phase masks of exp(iφm) in its respective Fourier domain O(kf). Through the forward and backward Fourier operation, a set of modulated images are averaged,

u(r)=1m∑m|um(r)/u¯m(r)|,

where

4um(r)=Fr−1[Fr(o(r))⋅exp(iφr,m)],

where u¯m(r) is the mean value of um(r), F denotes the 1-D Fourier transform, φ is normal random numbers with a mean of 0. RPM algorithm has two tuning parameters: the number of phase masks m, and the standard deviation of normal random numbers σ.

To quantitatively evaluate the performance of Deconv-RPM, we conducted 1-D simulations in line with the processing steps depicted in Fig. 1. FWHM and signal-to-noise ratio (SNR) are used as key metrics to compare the deblurring capability and signal-to-noise level between the unprocessed and the Deconv-RPM images. The SNR is defined as SNR = 10log₁₀(MAX²(u_img)/Var(n))²⁶, where MAX(u_img) represents the maximum pixel intensity of u_img, Var(n) denotes the variance of u_img in a background noise region n. The second row in Fig. 2(a) shows the noisy image u_img, generated by convolving ideal points (shown in the top row of Fig. 2(a)) with a complex Gaussian PSF (transverse) and then adding Gaussian white noise denoted as n. The FWHM and the SNR of intensity profiles were measured. The variance of the normally distributed noise, Var(n), was determined by a reasonable SNR value of 25 dB, approximated from the SNR of a conventional OCT image (refer to Fig. 3(a)). In estimating the real-valued PSF, we established a sub-optimal w, which is set to be 1.1 times the optimal w_opt. The tuning parameters m and σ were set to 10000 and 0.7, respectively. Upon examination of the unregularized deconvolved image (third row of Fig. 2(a)), we noted the presence of apparent ringing artifacts. In contrast, Deconv-RPM offers artifact-free images and demonstrates a 2.5-fold reduction in FWHM (evident in the bottom row of Fig. 2(a) and Fig. 2(b)). The process of image restoration by Deconv-RPM generally takes approximately 70 seconds for a matrix of 1000 pixels × 100 pixels. This process involves 10 iterations in the first step and 10000 resamplings in the second, executed on a 3.4 GHz Intel Core i7-3770 processor with 8 GB memory. Each resampling operation is performed independently, suggesting the potential for significantly faster computation through the use of GPU-based parallel computing.

We then evaluated how the parameters w, σ and m in the Deconv-RPM method contribute to FWHM reduction and SNR improvements. Fig. 2(c) illustrates a reduction in FWHM when w deviates from the optimal w_opt. When we varied the standard deviation σ, it was found that an optimal σ exists for both FWHM reduction and SNR enhancement, provided that other parameters are fixed (Fig. 2(d)). Fig. 2(e) demonstrates the algorithm’s performance reaching a limit as m increases. Under different noise levels in the unprocessed image (Fig. 2(f)), Deconv-RPM proved robust in FWHM reduction, with a SNR gain exceeding 10 dB in low SNR scenarios (unprocessed images, SNR below 32.5 dB).

The deblurring capability of the Deconv-RPM method was demonstrated through imaging of suspended polystyrene calibration spheres (80177, Fluka, diameter 2 µm) in water. A detailed description of the spectral domain μ-OCT system construction is presented in our earlier work^{27, 28}. A 20× objective lens (M Plan Apo NIR 20X, Mitutoyo, Takatsu-ku, Kawasaki, JP) was employed in the sample arm. The spatial resolution of μ-OCT in air was recorded as 2.40 μm × 1.38 μm (x, z). We observed high reflectivity signals, which originate from alterations in the refractive index (RI) at the peripheries of the spheres and the water interface. Particularly, these signals are manifested as paired reflections in the OCT images, each corresponding to the top and bottom surfaces of the microsphere, as illustrated in Fig. 3 inset. In out-of-focus regions, these paired signals widen due to the broadened transverse PSF.

Fig. 3(a, b), and Fig. 3(c) show logarithmic images following a two-step processing. Ringing artifacts, an expected outcome of the Richardson-Lucy deconvolution (Fig. 3(b)), were effectively mitigated by numerical random phase masks (Fig. 3(c)). The Gaussian kernels, derived from the theoretical physical parameters (the coherence length of the light source and the beam radius, respectively) and employed for deconvolution, were constant and non-adaptive. It's noteworthy that the numerical random phase masks, while suppressing ringing artifacts, also reduce some weak signals. We evaluated the system's performance based on the FWHM, under the assumption of a refractive index of 1.58 at 850 nm. When focused on the confocal region, Deconv-RPM achieved a 2.8 ± 0.2 times finer lateral FWHM and a 2.9 ± 0.5 times finer axial FWHM in comparison to conventional OCT. Two selected regions from Fig. 3(a) and Fig. 3(c) were magnified for detailed comparison of deblurring capabilities before and after Deconv-RPM processing. Two closely positioned spheres with four reflective signals, blurred in the conventional image, were clearly distinguishable in the Deconv-RPM image (Fig. 3(d)). Even at the depth range limit, Deconv-RPM's ability to distinguish reflective signals from merged ones highlights the improvement in image quality. The corresponding line profiles across reflective signals before and after Deconv-RPM processing are shown in Fig. 3(f) and Fig. 3(g).

Tissue experiments were conducted under the approval of the Institutional Review Board (IRB) of Nanyang Technological University (IRB-2016-10-015 and IRB-2019-05-050). The selection of parameters in the following tissue experiments was guided by the quantitative evaluation metrics derived from our simulations and phantom studies using microspheres. The theoretical FWHM w was set as the theoretical estimate of spatial resolution, and the number of iterations t for Richardson-Lucy deconvolution ranged between 4–10. Parameters σ and m were empirically selected to ensure artifact suppression. The pixel size in both spatial and frequency domains was set to fulfill the Nyquist sampling criterion.

To demonstrate the applicability of Deconv-RPM to static scatterers, we recorded the B-scans of swine floor of mouth exvivo using μ-OCT. In cross-sectional µ-OCT images, we usually identified paired high-scattering signals at the nucleocytoplasmic boundary each sandwiching a low-scattering nuclear core in the middle and upper epithelial layers. These corresponded to the Periodic acid-Schiff-Diastase (PAS-D) positive cells with flattened nuclei, as indicated in the matched histology. Speckle noise is inherent in coherent imaging²⁹, which reduces contrast and makes it difficult to resolve boundaries between highly scattered structures³⁰ in the tissue. For speckle noise reduction, we applied standard B-scan averaging (also known as spatial compounding) to both unprocessed OCT images and Deconv-RPM images. The principle underlying this operation is that the reduction in speckle noise is inversely proportional to the square root of the number of compounded frames (in this case, N = 20). However, there is an inherent trade-off between the speckle reduction and the resolution improvement: averaging B-scan comes at the expense of image resolution. This is primarily because the total field of view or angular aperture of the detector must be split to perform the averaging, inherently resulting in some loss of resolution³¹. Fig. 4 features cross-sectional imaging of swine floor of mouth: Fig. 4(a) is the standard μ-OCT, Fig. 4(b) represents the Deconv-RPM, and Fig. 4(c) is the deconvolution-only image. Magnified sections from Fig. 4(a) and Fig. 4(b) are detailed in Fig. 4(d) and Fig. 4(e), with nucleocytoplasmic boundaries marked by pink arrowheads. The histology reference from the swine mouth floor is provided in Fig. 4(f), and line profiles from Fig. 4(d) and Fig. 4(e) are shown in Fig. 4(g) and Fig. 4(h). Significantly, the Deconv-RPM depicted in Fig. 4(b) greatly improves the visibility of cellular structures. Nevertheless, because of the magnitude normalization implemented during the summation process in Step 2, the layered structures evident in Fig. 4(a) are less prominent in Fig. 4(b), a consequence of the self-averaging effect.

Dynamic imaging of cellular structures presents particular challenges due to the uncertainty of the physical model for digital correction and the abrupt changes in noise levels. To demonstrate the efficacy of our Deconv-RPM method in visualizing dynamic scatterers (Fig. 5 and Fig. S3), we performed labial mucosa imaging of a human volunteer using μ-OCT (a 126 kHz A-line rate allows for a B-scan rate of 512 lines at 246.1 frames per second). We consecutively scanned the same location on the mucosa multiple times to ensure the visibility of moving structures, such as red blood cells, in each B-scan. Fig. 5 presents blood flow in a capillary loop containing a mass of blood cells, where the flow direction is marked by red arrow. It should be noted that a single bright spot in the image does not necessarily represent an individual cell. Instead, it denotes a reflective surface. Fig. 5(c) was created by manually segmenting the capillary region in Fig. 5(a) and then applying Deconv-RPM to Fig. 5(b). To offer a detailed comparison between images processed with and without Deconv-RPM, we provided a 3× close-up view of the marked regions in Fig. 5(b) and Fig. 5(c) (Fig. 5(b') and Fig. 5(c')). This effectively highlights the real-time characteristics of blood cells. We ensured the registration of the capillary region by employing a single-step discrete Fourier transform (DFT) algorithm to align the B-scans³². The supplementary Movie S1 provides detailed x-z views of moving structures within invivo human labial mucosa as visualized by our Deconv-RPM method. Figure S3 presents an invivo imaging of zebrafish larvae captured through μ-OCT. An individual RBC can be seen in the tail artery, as highlighted in the orange box and inset. After applying the Deconv-RPM, the individual RBC becomes clearer and is easily distinguishable. The results presented above demonstrate that our proposed method can noninvasively enhance image quality for moving structures in real-time using μ-OCT.

To further understand the modulation of random phase masks, we conducted a series of 1-D simulations. Our initial assessment focused on the impact of normal distribution mean values on Deconv-RPM image quality. Following this, we investigated the effects of various distributions, such as uniform, exponential, and Poisson, on the quality of Deconv-RPM images. For consistent findings, each configuration was run 100 times independently. Each simulated matrix had dimensions of 1000 pixels × 1 pixel. The iteration number in the deconvolution process was 10, with the sub-optimal w set at 1.1 times the optimal w_opt. We fixed the resampling number m at 10000. This approach yielded a rich dataset from both original unprocessed and Deconv-RPM images, focusing primarily on SNR and FWHM metrics.

For the initial test, simulations were conducted within the range [−1, −0.5, 0, 0.5, 1] × π. Importantly, the standard deviation of the normal random numbers σ was fixed at 0.7. Fig. 6(a) illustrates that as the mean value of the phase distribution fluctuates, there is no significant change in the FWHM. Similarly, Fig. 6(b) indicates stable SNR values across different phase mask configurations. We used a one-way Analysis of Variance (ANOVA) test to determine differences among the five groups. The obtained p-values were 0.13 for the SNR of the original images, 0.25 for their FWHM; 0.81 for the SNR of Deconv-RPM processed images, and 0.48 for their FWHM. Given that all the p-values exceeded the standard significance level of 0.05, we could not reject the null hypothesis. This suggests that our analysis didn't identify any statistically significant differences among the group means for the metrics mentioned.

For the second test, we investigated various random phase distributions of phase masks. In the Gaussian distribution, we set the phase range as [−A, A]×π, with A as the variable. In the exponential distribution, we evaluated RPM modulation over varying phase means. In the Poisson distribution, changes in the rate parameter λ were examined. As shown in Fig. 7, the SNR (blue solid line) and FWHM (red solid line) of the original image remained stable, acting as reference values. The RPM-modulated values (dashed lines) showed deviations from these references based on the chosen distribution and its parameters. In Fig. 7(a) (uniform distribution), there was a peak in RPM modulated SNR and FWHM around the phase range of [−0.9, 0.9]×π, followed by a decrease. In the exponential distribution (Fig. 7(b)), the SNR increased with the phase mean, stabilizing between 1.7π and 3.7π. The RPM modulated FWHM became consistent after a phase mean of 2.5π. In the Poisson distribution (Fig. 7(c)), variations in both SNR and FWHM modulated by the RPM were observed, with λ = 0.95π identified as an optimal modulation point. In essence, each distribution presents unique modulation characteristics. Across all distributions, the relationship between SNR and FWHM remained consistent, suggesting a direct link between image clarity and resolution.

We attempt to provide a theoretical explanation and assume the effective signal to be a Gaussian function. There are two types of original signals under consideration: white noise and Gaussian signals. Furthermore, it's assumed that random phase masks follow a normal distribution.

Case 1: White noise signal f(x)

The procedure involves convolving the input signal f(x) with a random phase mask h(x). Formally, this can be represented as:

5c(x)=f(x)⊗h(x).

The autocorrelation function of h(x) is given as:

6Rh(τ)=1Nδ(τ).

Assuming that f(x) is white noise, its autocorrelation function R_f(τ) is a Dirac delta function:

7Rf(τ)=δ(τ).

Given that the autocorrelation function of the convolution c(x) is given by the convolution of the autocorrelation functions of f(x) and h(x), we have:

8Rc(τ)=Rf(τ)⊗Rh(τ)=δ(τ)⊗1Nδ(τ)=1Nδ(τ).

Case 2: Gaussian signal f(x)

Assuming f(x) is a Gaussian function represented by:

9f(x)=e−x22σ2.

To determine the autocorrelation function Rf(τ) of the Gaussian function f(x), we integrate the product of f(x) and f(x+τ) over all x:

10Rf(τ)=∫−∞∞e−x22σ2e−(x+τ)22σ2dx.

This evaluates to:

11Rf(τ)=πσe−τ24σ2.

Following the approach from Case 1, the autocorrelation function Rc(τ) of the convolution c(x) is given by:

12Rc(τ)=Rf(τ)⊗Rh(τ)=(πσe−τ24σ2)⊗1Nδ(τ)=πσNe−τ24σ2.

In summary, for white noise signals (Case 1), the autocorrelation function of c(x) is Rc(τ)=1Nδ(τ). The modulation with RPM results in an autocorrelation function that is entirely uncorrelated. This suggests that the noise remains at the level of a single image, even after superposition. For Gaussian signals (Case 2), the autocorrelation function of c(x) is Rc(τ)=πσNe−τ24σ2. The modulation with RPM leads to a correlated autocorrelation function. This indicates an amplification of signal strength when multiple images are superimposed. We could imagine that c(x) acts as an encoding process in an optical imaging encryption. For white noise, this encoding is stationary³³.

Applying deconvolution algorithms to noisy data, like that in OCT images, often poses challenges. While our approach has been successful in identifying cellular and sub-cellular features, there is some opportunity for further improvement. Some of the weak signals along with the artifacts generated by deconvolution could be also suppressed in step 2 of our approach, when they share similar SNR. This results in the unintentional loss of weak signals in our processed images. Addressing this overlap is essential, and several strategies can help in this differentiation. Analyzing the consistent behavior of pixels over time, known as temporal analysis, can help identify weak signals that stand out from the random fluctuations of ringing artifacts. Additionally, by examining the spatial patterns in the data, we can leverage the inherent structure of weak signals, which is typically absent in ringing artifacts. Analysis in frequency domain could also allow us to separate the unique components of weak signals from those of ringing artifacts. The integration of recent advanced deep learning algorithms offers another avenue, as these can be trained to discern even subtle differences between artifacts and weak signals.

Generally, when utilizing a high numerical aperture in the sample arm, OCT systems face limitations in depth-of-focus. Deconv-RPM is compatible with numerical refocusing methods such as interferometric synthetic aperture microscopy (ISAM)⁸, aperture synthesis³⁴ and subaperture correlation³⁵. In practice, scatterers outside the confocal region could be refocused before the Deconv-RPM processing. In strict terms, Deconv-RPM should depend on phase stability. Since the unstable phase fronts among axial scans do not greatly modify the backscattered intensity, Deconv-RPM is also robust to these errors as Richardson-Lucy deconvolution¹⁴.

The standard deconvolution, in concert with the numerical random phase masks, permitted artifact-free imaging with an approximate 2.5-fold FWHM reduction, which significantly exceeds the performance of these two methods separately. In conclusion, this method can be directly applied to conventional OCT images without any system modifications. It opens up new avenues for conducting deblurring imaging with potentially significant noise reduction, which could facilitate the identification of specific sub-cellular structures in turbid tissues. On a broader scale, the Deconv-RPM method can be applied to any other coherent or incoherent imaging modality suffering from severe noise artifacts, thereby yielding enhanced image quality.

This research is supported by the Guangdong Natural Science Fund General Program (2023A1515011289), Singapore Ministry of Health's National Medical Research Council under its Open Fund Individual Research Grant (MOH-OFIRG19may-0009), Ministry of Education Singapore under its Academic Research Fund Tier 1 (RG35/22) and its Academic Research Funding Tier 2 (MOE-T2EP30120-0001), and China-Singapore International Joint Research Institute (203-A022001).

The authors declare no competing financial interests.

Supplementary information for this paper is available athttps://doi.org/10.29026/oes.2024.230020