Non-Local Mean Denoising Algorithm Based on Fractional Compact Finite Difference Scheme Effectively Reduces Speckle Noise in Optical Coherence Tomography Images

Optical coherence tomography (OCT) is used in various fields such, as medical diagnosis and material inspection, as a non-invasive and high-resolution optical imaging modality. However, an OCT image is damaged by speckle noise during its generation, thus reducing the image quality. To address this problem, a non-local means (NLM) algorithm based on the fractional compact finite difference scheme (FCFDS) is proposed to remove the speckle noise in OCT images. FCFDS uses more local pixel information when compared to integer-order difference operators. The FCFDS operator is introduced into the NLM algorithm to construct a high-precision weight calculation so that the proposed algorithm can effectively reduce the speckle noise in the OCT images. Experiments on simulations and real OCT images show that the proposed method is comparable to other state-of-the-art despeckling methods and can substantially reduce noise and preserve image details such as edges and structures. Speckle noise removal can further promote the application of the proposed algorithm in medical diagnosis and industrial detection, as it has key research value.

Optical coherence tomography (OCT) is a non-invasive, high-precision imaging modality that is used in various fields, such as dermatology and ophthalmology [1]. However, as the OCT system uses the principle of light interference, the generated image is severely affected by speckle noise. This leads to a significantly low-quality image [2]. The layered structure and fine features of OCT images are blurred due to their low quality, which affects the accuracy of the diagnosis by the doctor [3]. An effective OCT image denoising algorithm can improve the image quality, thereby enhancing the quality of the diagnosis.
Recent researches show that either direct averaging [4] or hardware-based phase modulation [5] then averaging multiple OCT images can significantly suppress the speckle noise and enhance the OCT image quality with the cost of longer time data acquisition. However, algorithm-based speckle denoising methods are still in high demand. In addition, many despeckling algorithms have been proposed for OCT images, which are generally divided into two categories: Transform domain method and spatial domain method. Transform domain algorithms such as the wavelet shrinkage technique [6,7] and the curvelet shrinkage technique [8] usually reduce the speckle noise in OCT images by manipulating the size of the correlation coefficients. However, it is difficult to distinguish fine details from the speckle noise using the transform domain algorithms. Some artifacts may also be produced during the denoising process. Spatial domain algorithms are usually divided into local and non-local filters, such as partial differential equation (PDE) methods [9,10], nonlocal mean (NLM) methods [11][12][13], sparse representation (SR) methods [14][15][16][17], low-rank approximation methods [3,[18][19][20], multi-frame dependent algorithms [21][22][23][24] and deep learning methods [25][26][27][28]. The deep learning method is a relatively advanced algorithm that requires a large amount of training data. The training of the deep learning models requires high-performance computing and is time-consuming. In addition, the estimated image obtained by the PDE algorithm is under-smooth, while the estimated image obtained by the SR, LRA, or multi-frame dependent algorithms are over-smooth. Hence, the NLM algorithm has attracted the attention of many scholars.
The NLM algorithm reconstructs image pixels by the weighted averaging of pixels in a predefined search window, where non-local weights are determined based on the similarity between different image patches rather than individual pixels. The NLM method, originally designed for Gaussian noise removal, has lately been applied to remove speckle noise from ultrasound images [29,30] and synthetic aperture radar images [31]. Furthermore, by designing a new weight calculation operator, some NLM-based filter filters are proposed for OCT image denoising. A double-precision NLM filter with a Gaussian anisotropic kernel function, instead of the traditional uniform kernel function, was previously proposed to effectively remove speckle noise [11]. Additionally, a two-stage NLM method using the uncorrupted probability of each pixel was introduced to suppress noise [12]. In literature [32], they have presented a novel despeckling technique for OCT, TNode, which is based on the non-local means algorithm with proper speckle statistics. This algorithm could remove the noise in the background, but it was not sufficient to remove the noise in a high-intensity region. In addition, an NLM method based on guided filtering was also proposed for the despeckling of the OCT images [13]. This algorithm could remove the noise in the background, but it was not sufficient to remove the noise in a high-intensity region [33]. An NLM method based on guided filtering was also proposed for the despeckling of the OCT images. However, these NLM-based methods could only marginally improve the performance of the despeckling system. The disadvantage of NLM methods is that they cannot accurately determine the weights using only the neighboring individual pixels. Hence, the structural similarity of noisy image blocks cannot be effectively represented. Based on existing research, it can be understood that further characterizing the topological structures as image edges and textures by substantially reducing noise and improving the robustness of the algorithm is a very critical scientific issue. To address this problem, we propose an improved NLM algorithm, which uses a fractional compact finite difference scheme (FCFDS) to accurately calculate the image similarity weight. The proposed algorithm can provide a better filtering effect, as it not only fully removes speckle noise but also effectively retains the topological structure of the image. The enhanced filtering effect can further promote the application of the OCT image data in medical diagnosis, industrial inspection, cultural relic restoration, and other fields. The improved despeckling algorithm has significant research value. The results of this study can provide important technical support for the early diagnosis and treatment of diseases. The primary aims of our research are: (1) A template convolution operator based on the FCFDS is proposed herein. The FCFDS is a high-precision format, an extension of the compact difference format for integerorder equations. Generally, a format with higher precision implies the use of more points and a greater number of calculations. However, the FCFDS can achieve high accuracy with fewer points. Hence, using the FCFDS can reduce the computational complexity of the algorithm while maintaining its effectiveness.
(2) A bidirectional, high-precision image similarity weighting strategy based on the FCFDS is used. The information regarding the main structure of the OCT image is primarily concentrated in the horizontal and vertical directions. Therefore, to better preserve its structural information, we propose a bidirectional non-local weighting operator. In addition, we propose to use the high-precision FCFDS operator to construct the variance parameter in the vertical direction to improve the preservation of the detailed texture of the image.
(3) An OCT image denoising algorithm based on the FCFDS and NLM algorithm is proposed to remove noise while effectively preserving the structural information of the image.
The rest of the paper is organized as follows. Section 1 provides the detailed process of the proposed method. In Section 2, the simulation and real OCT images are studied, and the experimental results obtained are discussed. Section 3 concludes the paper.

Proposed Method
The main advantage of fractional differentiation over integer order differentiation is that fractional-order differential operators are non-local and provide an excellent tool for describing the memory and genetic characteristics of image processing. However, in the classic integer order model, the influence of fractional order differentiation is ignored [34]. In addition, the higher-order calculus derivatives and the generalized differential quadrature method can obtain higher-precision orders in the numerical solution of partial differential equations. These higher-precision orders are more conducive for applications in practical engineering calculation problems. However, the construction of higher-order derivatives or generalized differentials requires more mesh points, increasing the complexity of the algorithm and its computation Fractional differential operators have also achieved good filtering effects in image denoising, but they have fewer applications in OCT image denoising. In addition, the development of the FCFDS increases prospects for applications of the fractional derivative theory and improves the accuracy of calculations. This supports fractional calculus in the preservation of non-local image block detail information. The following two lemmas are used in the difference scheme to derive the FCFDS.
where α ∈ (0, 1), β, γ, and k are constants uniformly for x ∈ R as γ −→ 0, and o(β) represents the higher-order infinitesimal of β. −∞ D α x f (x) is the Riemann-Liouville fractional derivative, which can be expressed as: where ε is the integral variable, and n is the smallest positive integer greater than α.
x f (x) and its Fourier transform belong to L 1 (R). Defining the weighted and shifted Grünwald difference operator by , and γ, p, q are all integers. Then, we have uniformly for x ∈ R as γ −→ 0. Because of the symmetry of γ, p and q, it may be assumed that δ > p > q. By choosing γ = 0, p = −1 and q = −2, we get According to Lemmas 1 and 2, if } is the grid function space, then for any grid function f ∈ V β , the FCFDS operator is defined as: where n = 2, 3, · · · , N.
From Equation (6), we can see that the FCFDS operator can be simplified to approximate expressions for multiplication and addition. When β = 1t for a two-dimensional image f (x, y), the corresponding expressions in the horizontal x and vertical directions y directions are shown in Equations (7) and (8).
Since the structural information of the OCT images is concentrated in the horizontal and vertical directions, fractional differential masks are constructed in the directions of the positive x-coordinate and positive y-coordinate, respectively, as shown in Tables 1 and 2. Note that c k , k = 1, 2, · · · , n is the mask coefficient of the pixel of interest. When k ← N = 2s − 2, the fractional differential mask (2s − 1) × (2s − 1) can be implemented. To make sure fractional differential mask has the certain center, in general, an even number is selected as N. It can be seen from Equations (7) and (8) that the mask coefficient c k , k = 1, 2, · · · , n can be expressed as: The construction of the traditional mean operator leads to the loss of image details. To better preserve the details of the image, based on the FCFDS operator, we calculate a new similarity weighted kernel function. The new weight value can be defined as where v(•) is a pixel value at the position • from a noisy image v; w(x, t) is used to evaluate the similarity between the pixel intensity of two image blocks, where x and t represent two square image blocks of size k × k. As a filtering parameter, h controls the degree of decay of the exponential function. In addition, ∆ can be represented as (∆ x , ∆ y ), where (∆ x , ∆ y ) ∈ ∆ and ∆ ∈ Θ; Θ can represent the image block region; (∆, x) represents the kernel function. Unlike the traditional convolutional Gaussian kernel function, the new kernel function (∆, x) is constructed based on the FCFDS operator and contains more non-local information. (∆, x) can be written as: where σ x ∈ [−k, k] and σ y ∈ [−k, k]; φ x and φ y (∆, x) represent the variance parameters of the image block in the horizontal and vertical directions, respectively. Notably, the φ x is a constant, and φ y (∆, x) is determined by the FCFDS operator, and its expression is: where η is a damping coefficient, which determines the sensitivity of the φ y (∆, x) value depends on ∆ α x and ∆ α y .
In this study, we use the weight function constructed by Equation (10) to obtain the filtered result image by weighted average, which is defined as where F(x) is a set of neighborhood pixels around x, which are referred to as searching window.
To summarize, the pseudo-code of the proposed method is described in the Algorithm 1.

Algorithm 1
Proposed method for OCT image denoising.

Experimental Results
In order to verify the effectiveness of the algorithm, we perform noise removal experiments on real OCT retinal images, and some representative test images are shown in Figure 1. Where (a)-(c) and (d)-(f) in Figure 1 were provided online by the authors of Reference [15] and [38], respectively. Note that the DUKE dataset in Reference [15] and the PKU37 dataset in Reference [38] provide clean label images, that is, the average image of multiple adjacent images. In addition, for performing visual inspection and a quantitative comparison, the proposed algorithm was compared and analyzed against the outputs obtained after using NCDF [9], PNLM [12], WGLRR [19] and PBFD [3] four more advanced denoising algorithms. The quantitative indicators used include Peak Signal-to-Noise Ratio (PSNR) [39] and Structural Similarity Index Measurement (SSIM) [39]. PSNR is mainly used to measure the sharpness, contrast and clarity of images. The larger their values, the better the structure of the image such as edges and textures is preserved. SSIM is mainly used to measure the smoothness of the image. The larger their value, the greater the smoothness of the image. All the algorithms are implemented using MATLAB programming on a laptop with the following hardware specifications: 3 GHz CPU and 16 GB RAM.

Visual Analysis of Real OCT Image Despeckling
The speckle noise removal experiment was carried out on the real OCT retinal test image, and the denoising effects of NCDF, PNLM, WGLRR and PBFD algorithms were compared. Some visual quality results after denoising are shown in Figures 2 and 3. From Figures 2 and 3, it can be seen that the denoising result of NCDF algorithm is compromised. Although the NCDF algorithm reduces some speckle noise on OCT image, its structure information is still blurred. In addition, the PNLM algorithms can fully remove the noise in the background area, but the high-intensity structure area still contains some noise and the layer structure information is not clear. Furthermore, one finds that the WGLRR algorithm and the PBFD algorithm can fully remove speckle noise on OCT images, but over-smoothing may lead to layer structure confusion in the edge region of the filtered image. Finally, one sees that the filtering effect of the proposed algorithm is similar to that of the WGLRR and PBFD algorithms. It can not only reduce the speckle noise in the image, but also preserve the strong and weak edge structure information and detail information of the image. The estimated image obtained by the proposed method can well preserve the edges and features in the OCT image, which is very useful for other medical image analysis, such as edge detection, image layer segmentation and clinical diagnosis.

Quantitative Analysis of Real OCT Image Despeckling
In order to more clearly show the robustness and superiority of the proposed algorithm, in addition to the previous visual effect comparison analysis, we also quantitatively analyze the filtered images of the proposed algorithm and the four comparison algorithms, and the results are shown in Tables 3 and 4. Note that since the clean label image given in the dataset is obtained by weighted average of multiple frames, there is a certain error in the calculation of PSNR and SSIM values, but we believe that the error is acceptable.
From Table 3, one finds that the PSNR values of the proposed algorithm for DUKE and PKU37 test images are greater than the PSNR values of the NCDF and PNLM algorithms, so the smoothing effect of the proposed algorithm is better than the NCDF and PNLM algorithms. In addition, one can also see that although the PSNR mean of the proposed algorithm is lower than the PSNR mean of the WGLRR and PBFD algorithms, he difference between their values is small and the PSNR value of the proposed algorithm for some test images in the DUKE dataset is better than these two algorithms. Therefore, the smoothing effect of the proposed algorithm for DUKE test images is similar to that of WGLRR and PBFD algorithms. From Table 4, it can be seen that the SSIM values of the proposed algorithm for DUKE test images are similar to the results of PSNR. The SSIM values of the filtered images of the proposed algorithm are greater than those of NCDF and PNLM algorithms. The SSIM values of individual images have achieved the best results, and the SSIM values of some filtered images are greater than those of the PBFD algorithm. Therefore, for the Duke dataset, the ability of the proposed algorithm to preserve the structure is better than that of NCDF and PNLM algorithms, which is comparable to WGLRR and PBFD algorithms. In addition, for the PKU37 dataset, although the proposed algorithm does not obtain the highest SSIM values, their values are not the worst, that is, the proposed algorithm can not only better reduce speckle noise but also better preserve the structural information of the image.  In summary, the proposed methods for different test images are robust and have better denoising effect. Such an output image can effectively promote its subsequent application in scientific research.

Analysis of Main Contributions
The main contribution of this paper is to propose a new kernel function based on FCFDS format, and use the new kernel function to construct high-precision weighted operator. In order to better demonstrate the effectiveness of the proposed new kernel function, we used the same algorithm and parameters to conduct denoising experiments with the traditional kernel function and the new kernel function respectively. The experimental results are shown in Figures 4 and 5.  It can be seen from Figure 4 that the new kernel function can greatly improve the filtering results of OCT images. The filtering result of the NLM algorithm using the traditional kernel function is under-smooth. The proposed algorithm using the new kernel function can effectively reduce the speckle noise of the OCT image and preserve its structural information. In addition, from Figure 5, one sees that the PSNR and SSIM results obtained by the algorithm using the new kernel function are much larger than the NLM algorithm using the traditional kernel function. Therefore, the validity of the new kernel function is demonstrated again from the perspective of quantitative index analysis.

Visual Analysis on OCT Data
The OCT estimated image after speckle removal is an important tool for post-processing such as edge detection and segmentation. In order to further evaluate the effectiveness of the proposed algorithm, we use Canny edge detection algorithm [40] to test the edge of the denoised image. The results of image edge detection after denoising by each comparison method are shown in Figure 6.
From Figure 6, one finds that all comparison algorithms except the NCDF algorithm can better detect strong edge information on the image. However, the filtering result of WGLRR algorithm is too smooth, weak edges and details are difficult to detect. In addition, similar to clean label images, one also observes that the PNLM algorithm, FBFP algorithm and the proposed algorithm can all detect strong and weak edges and details, but it can be seen from Figure 6 that the noise in the background region of the PNLM algorithm and FBFP algorithm is also detected. In short, through edge detection, it can be found that the proposed algorithm has well removed the noise of the image, and clearly displayed the layer structure of the image, and the weak boundary layer has been well preserved.

Conclusions
We proposed an NLM algorithm based on the FCFDS to reduce the speckle noise in OCT images. The proposed algorithm primarily utilizes the non-local self-similarity, statistical properties, and FCFDS operator theory of the image itself. To examine the validity of the approach applied in this study, the experimental results were compared with those published in the previous research, and a good agreement was observed between them. The experiments results revealed that: The bidirectional high-precision image similarity weighting strategy based on the FCFDS is more effective than the traditional isotropic Gaussian similarity weighting strategy for denoising the OCT images.
Compared to some commonly used methods to remove speckle noise, our proposed method obtains the better results in both visual inspection and quantitative index analysis.
Therefore, the denoising algorithm proposed herein is expected to be applied in a wide range of biomedical imaging applications and to promote scientific research related to its subsequent processing.