Adaptive anisotropic diffusion for noise reduction of phase images in Fourier domain Doppler optical coherence tomography

: Phase image in Fourier domain Doppler optical coherence tomography offers additional flow information of investigated samples, which provides valuable evidence towards accurate medical diagnosis. High quality phase images are thus desirable. We propose a noise reduction method for phase images by combining a synthetic noise estimation criteria based on local noise estimator (LNE) and distance median value (DMV) with anisotropic diffusion model. By identifying noise and signal pixels accurately and diffusing them with different coefficients respectively and adaptive iteration steps, we demonstrated the effectiveness of our proposed method in both phantom and mouse artery images. Comparison with other methods such as filtering method (mean, median filtering), wavelet method, probabilistic method and partial differential equation based methods in terms of peak signal-to-noise ratio (PSNR), equivalent number of looks (ENL) and contrast-to-noise ratio (CNR) showed the advantages of our method in reserving image energy and removing noise.


Introduction
Since the invention of optical coherence tomography (OCT) in the 1990s [1], it has gradually become a powerful imaging modality in the field of medical diagnosis [2][3][4]. Combined with Doppler technique, simultaneous structure and flow imaging provided by phase resolved Doppler optical coherence tomography (PRDOCT) expanded the applications to flow research, such as blood vessel [5], retinal blood flow [6,7], cerebral blood flow [8], and blood velocity and flow from tissues [9]. High quality and clean phase images are the basis for not only flow speed extraction but also various phase based postprocessing techniques such as phase unwrapping, image segmentation and three-dimensional flow field simulation. Unfortunately, phases images from OCT system usually suffer from serious noises, which makes effective noise reduction methods for phase images highly desirable. However, a lot of attention has been paid to structural images of OCT. There is a lack of methods available specially designed for phase images of OCT.
OCT images mainly contain speckle noise and random noise. In early OCT noise studies, lots of attention were paid to the noise property and formation mechanism in OCT system, especially in speckle noise. J. M. Schmitt et al. [10] explained that speckle has double property: signal carrier and noise factor. Due to the coherent imaging property of OCT, speckle noise becomes inevitable. With the increased understanding of noise property, various algorithms have been developed for noise reduction.
These algorithms can be classified into four groups: filtering method, wavelet method, probabilistic method and partial differential equations (PDEs) technique. Typical examples of the first group are mean filtering [11] and median filtering [12]. However, these methods cannot work well for low SNR images. The second group refers to wavelet method, which is usually combined with other methods. For instance, wavelet with Lee filtering and Wiener filtering were used to reduce speckle noise of the OCT images [13]. Wavelet denoising of multiframes have been utilized for denoising OCT data [14]. In the thrid group, Alexander Wong et.al utilized Bayesian estimation to reduce speckle noise [15]. The fourth group are based on PDEs, which have been widely used in image processing and computer vision. Back in 1990, Perona [21]. reduced speckle noise of OCT images by total variation regularization. In conclusion, diffusion model methods based on PDEs have attracted a lot of attention in the field of OCT.
However, most of methods mentioned above were applied to intensity images, while phase images are rarely processed. On the other hand, noise level of phase image is larger than intensity image. Object of this work is to develop an effective technique to remove noises in phase images. For the diffusion methods, diffusivity is highly dependent on image gradient. Due to the fact that noise has similar gradient thus diffusivity compared to signals at edges, original diffusion equation cannot handle noise very well. To address this issue, we propose an effective noise estimation criteria with two estimators to differentiate noise from signals more accurately. Different pixels based on the estimation criteria are assigned different diffusivity. What's more, 1/4 of the thermal diffusivity is used as the iteration step, which makes noise and signal have different iterative process for an image.
For the organization of this paper, Section 2 describes two noise estimators, noise estimation criteria and the adaptive anistropic diffusion model. Then PSNR, ENL and CNR are used to analyze experiment results and the diffusion process in Section 3, which includes the optimization process of our algothrim for phantom images and results from mouse artery images. Finally, we conclude in Section 4. Figure 1(a) and 1(b) are intensity image and phase image of a transparent plastic tube acquired from a 1300nm band OCT system, which has an axial resolution of 14 μm and imaging range of 6.7 mm. The OCT system was running at 70 fps with 1000 Ascans per frame. From the visual view, phase image presents many outliers in the signal and background region. To illustrate the characteristics of noise, two local regions are selected and enlarged in the right side of OCT images. From the comparision of intensity image and phase image in the box region, we can see that phase image has larger noise level than intensity image. For box 1, phase image has clear visualization of the random pixels with large fluctuation. In addition to random noises, box 2 of phase image has many speckle pixels with low amplitude fluctuation and this speckle property has reduced the smoothness of homogenous region. Thus we can see that random noise and speckle noise are main parts of the phase image noise. They share the same feature that central pixels have large difference from neighboring pixels. Therefore, this degree of difference can be used to estimate the noise probability of every pixel.

Local noise estimator(LNE)
To estimate noise, local noise estimator (LNE) [22] for a pxiel p is defined by the expression where ( ) 0 N q is 5 5 × neighbor pixels centered at pixel p excluding pixel p .
( ) p I q describing the relationship between the central pixel p and the 24-neighborhood pixels q can be calculated quantitively as where ( ) ( ) , u p u q are pixel values of , p q . And l T is the setting threshod. Equation (2) describes number of the pixel, which are similar to the central pixel. In reality, noise pixels are of larger gray value than those of its neighbor pixels, which means the isolated noise pixels have large LNE values and pixels uncorrupted by noise have small LNE values. In fact, LNE estimator judges the similarity between the central pixel and its surrounding pixels, thus it is a local estimator. Parameter l T is applied for image pixels, thus l T can be expressed as 2 r π , where r denotes the percentage of setting for all the pixels and 2π denotes range of phase value. When r is given, parameter l T is determined, thus LNE is a rigid estimator.

Distance median vaule (DMV)
To increase estimation accuracy, DMV estimator [23] is defined by where ( ) ( ) median I p denotes median value of 7 7 × neighbor pixels for a pixel p . Since the intensity of noise pixel is different from surrounding pixels, specifically the noise pixel is an outlier, DMV value of a noise pixel is larger than the the pixel uncorrupted by noise, thus DMV also becomes indicator for noise. From Eq. (3), we can see that DMV is also a local estimator. Due to the flexibility of ( ) median I , DMV is a flexible indicator compared to LNE estimator.
We determine the window size by a quick check of the LNE and DMV images with different window sizes first and then choose the one with sharper contrast with a prone towards small window size to reduce time cost.

Noise estimation criteria
In subsection 2.2 and 2.3, LNE and DMV estimator have been utlized to estimate noise. In order to unify these two estimations into diffusion coefficient, which is denoted as c , we propose a noise estimation criteria illustrated in Fig. 2. Here, 1 2 , T T are two thresholds corresponding to LNE and DMV for pixel identification. Based on the fact that LNE and DMV values of noise pixels are larger than signal pixels, we can assume that noise pixels correspond to . From the definition criteria, one can see that image pixels are classified into four groups and diffusion coefficient is also divided into four intervals. Signal pixel has small diffusivity and noise pixel has large diffusivity. Here, the setting interval point has a constraint 1 2 3 Computing process of the diffusion coefficient c is shown in Fig. 2. For a pixel, we need to determine its class first. Then, normalized weight 1 2 , ω ω of LNE and DMV for this pixel will be calculated depends on its category. After that, the total weight of these two weights denoted as ω is calculated. When the setting interval is [ From the calculation process for diffusion coefficient, we can see that every pixel is assigned a unique diffusion coefficient by determining the pixel class and the setting interval points. In reference [19], it has a similar computing process with an average weight between parameters. However these parameters are calculated with respect to all pixels, so signal pixel may have large diffusion coefficient. Thus, it is difficult to ensure a reasonable estimation of noise for this simple averaging operation. In Fig. 2, four types of pixels based on LNE and DMV promise that the diffusion coefficient divided into the reasonable interval, which means that our design of the diffusivity can achieve better performance. Values of these parameters ( ) , , , , , , x x x α ω ω ω will be explained in subsection 3.4.

Denoising model based on adaptive anisotropic diffusion
The general anisotropic diffusion equation is introduced to decribe the image diffusion process as follows [16]: where I ∇ denotes image gradient and ( ) c ⋅ denotes diffusion coefficient in Fig. 2. Then, a discretized approximation by the forward and backward differences shown as follows [24]: (6), we can see that noise pixel has strong diffusion action and signal pixel has weak diffusion action. Thus noise can be removed and signal will be kept.
It is worthy to note that many diffusion models adopt the constant step size [25-28] for each iteration or whole iterative process of the image. Here a better iteration step is propsoed in the Eq. (7).
where 1/4 is used to promise the convergence of the Eq. (5) [16]. Since diffusivity of noise pixel is large, the corresponding step is large and then this pixel can be removed quickly. Conversely, pixel uncorrupted by noise have a small diffusivity, the corresponding step is small and this pixel can be kept very well, thus every pixel has different step size, which depends on the diffusivity of this pixel. Therefore, Eq. (5) shows two adaptive property in the diffusion filtering process. The whole filtering process for phase image is shown in Fig. 3. For phase images with noise, when the parameter l T is given, we can calculate the LNE and DMV values using Eq.
when IE is less than or equal to tolerance ie T , the iterative process is stopped. In this paper, tolerance ie T is set as 0.003. Fig. 3. The implement filter process for a image.

Image quality evaluation
To compare results of different image processing methods, we select peak signal to noise ratio (PSNR) where , d o I I represent the processed image and the original image. h denotes pixel number of the image. Considering the fact that there is no original standard phase image, we use 20 frame images and the subpixel image registration algorithm [34,35] to make a reference image (Fig. 4(a)) for PSNR calculation.
ENL: ENL describes the smoothness of regions of interest (ROIs) corrupted by noise and the corresponding expression is u σ are pixel mean and standard deviation of background region.

Experimental results and discussions
In this section, our proposed method and other seven methods. PM diffusion [16], NCDF method [20], mean filter [11], median filter [12], wavelet method [12] and Bayesian method [15,36] are used to process the phase image with noise (shown in Fig. 1(b)). Then, we compare their results from the visual view and the index evaluation (PSNR, ENL, CNR and the execution time with 10 times averaged for various methods). For the experiment, image process was conducted through MATLAB R2013 in a computer with Intel i7 processer of 2.6 GHz and 8GB RAM.

Results of the visual comparison
First, we get processed images by our proposed method and other methods in Fig. 4. Figure  4(b) shows processed image after PM diffusion filtering. From the visual view, we can cearly see that it almost has no effect on noise. Figure 4(c), 4(d), 4(f), 4(g)) show the OCT images after NCDF, mean, wavelet, Bayesian filtering. They have poor ability in keeping image energy of ROIs because the maximum pixel value of the processed image is reduced to half of Fig. 4(a). From Fig. 4(e), we can see that median filtered image cannot remove noise completely although it can keep the maximum value of the pixel. However, our proposed method can overcome these problems. From Fig. 4(h), one can see that noise pixels are removed completely and center region obtains the best smoothness. More significantly, image energy is preserved to the maximum extent. Maximal phase value decreased from 3.14 to 2.91 for our method. Thus, the visual comparison clearly indicates advantages of our method in keeping the image energy and removing noise.

Results of the quantitative comparison
In order to be more persuasive, we use PSNR to describe the retention degree of image energy and LNE, CNR are used to describe the extent of noise smoothness. To calculate LNE, CNR values, eight regions of interest (labeled 2 to 9) and one background region (labeled 1) are selected and marked with black-line boxes in Fig. 5. The results are shown in Table 1. PM diffusion obtains the worst results (the smallest PSNR, LNE and CNR). The use of NCDF, mean, wavelet and Bayesian filtering result in some improvement in the image PSNR (15.26dB, 16.13dB, 15.22dB, 14.90dB and 15.29dB), LNE and CNR. The application of our filter results in a significant image quality improvement (PSNR is 17.43dB), which is also apparent by comparing the change between Fig. 5(h) and 5(a). The PSNR improvement with the our method is better than other methods. Besides, LNE and CNR also achieve the largest values, which means noise is removed well.   In Subsection 2.5, iteration error shown in Fig. 6(a) is introduced to show the change between the two adjacent images in the iterative process. From Fig. 6(a), one can see that the iterative error drops sharply at the beginning of iteration and then stablizes when the iteration is greater than 30. In fact, the iterative process can be divided into two parts. The first part mainly remove outliers, such as random noise. The second part corresponds to the process that uneven regions corrupted by speckle are smoothed.
Profiles of the red line marked in Fig. 6(b) are plotted to compare the smoothness of different methods. The curve using PM diffusion is basically coincident with the noise curve, which shows that this method cannot smooth noise. For NCDF, mean, Bayesian and wavelet methods, due to excessive filtering to image, their curves are not as high as our curve on the wrapped region of phase. Their curve is of larger fluctuation than ours. And for median, the curve fluctuation is larger than ours although meadian can keep large phase value.

Comparison of three classes diffusivities
In Subsection 2.4, we have explained the advantages of our method for calculating the diffusivity from the theoretical point of view. In order to show its advantages more directly, the normalized LNE and DMV are used separately as the diffusivity for comparison. In this comparsion, all conditions are the same as our method except for the diffusivity. The results are shown in Fig. 7. As expected, we can see that the maximum PSNR by using our diffusivity is larger than the other two ways, so these results again indicate that our method performs better. Smoothness of the Fig. 7(c) is better than Fig. 7(a) and Fig. 7(b) from the visual view.

Parameters optimization based on genetic algorithm
During the diffusion process, we need to set seven parameters ( ) To achieve the optimal results, genetic algorithm [37] is used to optimize these parameters and the optimization results are shown in Table 2. Figure 8 shows PSNR value of the optimization process. When the iteration number exceeds 50, PSNR value with small fluctuation is close to 18 dB.

Results of vessel phase images
In the above part, our method was applied to sample phase images and we can get the optimal result. In order to further perform the application of our method in medical field, we select three groups vascular images ( Fig. 9(a), Fig. 11(a), Fig. 13(a)) from the various positions of mouse blood vessels in OCT system reported in [38] and use seven methods to remove phase noise of these vessels. Figure 9, Fig. 11 and Fig. 13 are the processed results for visual comparison. Figure 10, and 12 are labeled with regions of interest and background for parameter comparison based on Fig. 9 and Fig. 11 respectivley. Here, please note that our method still use the optimization parameters of Table 2 to process the vessel phase images. Then, ENL and CNR are used to measure filtering results of different methods (Table 3, Table 4, Table 5). These results demonstrate that our method still gains attractive results on no matter which group of vascular images. Therefore, it is appropriate that the optimization parameters of sample phase image can be applied to these vascular images due to the similar noise property, although sample images and vascular images come from different OCT systems.  Based on Fig. 9, we seleted four regions of interest (Lable 2 to 5) and one background (Lable 1) shown in Fig. 10 to calculate CNR and ENL using different methods. Three of the four regions of interest are at the edge of sample and one is at the center of the sample.  Similarly, based on Fig. 11, we seleted four regions of interest (Lable 2 to 5) and one background (Lable 1) shown in Fig. 12 to calculate CNR and ENL using different methods.  Using the same approach as Fig. 10 and Fig. 12, we calculate the CNR and ENL based on Fig. 13 and results are shown in Table 5. For our method, time cost mainly includes three parts: the noise type estimation, anisotropic diffusion and parameters optimization. Parameters optimization is the most timeconsuming part with orders of magnitude higher than the other two parts. There are mainly two reasons. First, there are abundant of possible parameter sets to be evaluated. Second, the iteration number is intentionally fixed at large number for optimization for each parameter set. However, the optimization process needs to be done only once. Based on our evaluation, noise type estimation costs around 92.66% and anisotropic diffusion costs around 7.34% of the computation time for each iteration. Depending on the sample image, the total time cost varies from 17.54 seconds to 18.15 seconds.

Conclusion
In this paper, we present a noise estimation criteria by two estimators (LNE and DMV) for noise estimation in OCT phase images. Both theoretical explanation and experimental data show advantages of our noise estimation criteria. In addition, this estimation criteria is used to design a new the diffusivity and adaptive iteration step for the anisotropic diffusion, so the adaptive anistropic diffusion can effectively eliminate phase noise. It is an important topic to determine parameters 1  2  3  1  2  3 , , , , , , , r T T T x x x α for our method. Therefore, the genetic algorithm is used to optimize these parameters. PSNR, ENL and CNR are used to evaluate our method and other methods. Experimental results demonstrate that our method obtains the optimal results regardless of maintenance of the the maximum of the pixel phase value or the noise removal for sample phase images. More importantly, we apply the optimization parameters of sample phase image for vascular images and also achieve similar results, which indicates that we can use the parameters after the optimization of our method to other OCT phase images.
The major limitation of current method is relative high time cost. Nevertheless, this issue can be solved by utilizing parallel computing method such as general computation from graphics processing unit since the main time consumer, calculation of LNE, DMV and anisotropic diffusion can be implemented independently with each other at every pixel. This will be one of our future work to do.