Improved wavelet hierarchical threshold filter method for optical coherence tomography image de-noising

Jing Cao*, Pinghe Wang*, Bo Wu, Guohua Shi*‡**, Yan Zhang, Xiqi Li¶||, Yudong Zhang¶|| and Yong Liu* *School of Optoelectronic Information University of Electronic Science and Technology of China No. 4, Section 2, North Jianshe Road Chengdu 610054, P. R. China College of Optoelectronic Technology Chengdu University of Information Technology Chengdu 610225, P. R. China Suzhou Institute of Biomedical Engineering and Technology Chinese Academy of Sciences Suzhou 215163, P. R. China School of Electronic and Communication Engineering Guiyang University Guiyang 550005, P. R. China ¶Institute of Optics and Electronics Chinese Academy of Sciences Chengdu 610209, P. R. China ||Chinese Academy of Sciences The Key Laboratory on Adaptive Optics Chengdu 610209, P. R. China **ioe_eye@126.com


Introduction
Optical coherence tomography (OCT) is a promising optical imaging technique widely used in ophthalmology and cardiology, due to its advantages including non-invasive, high speed, high resolution and three-dimensional imaging. [1][2][3] The quality of images is vital for OCT. In reality, since OCT is based on low-coherence interferometry, [4][5][6] it inevitably brings speckle noise to OCT images. Speckle noise degrades OCT images' quality and makes the details of the image blurry. 7 For better understanding the speckle noise in OCT, many groups have done signi¯cant research on its statistics and properties. [8][9][10][11] Many methods have been proposed to attenuate the speckle noise in OCT images, these methods can be divided into two parts: hardware and software methods. Hardware methods are optical approaches that physically remove speckle noise, including frequency 12,13 and spatial compounding, [14][15][16][17][18] which can signi¯cantly improve signal-to-noise ratio of OCT images. But these methods would increase system complexity, meanwhile decrease imaging speed and spatial resolution. 12,16 Software approaches based on post-processing¯ltering techniques, including adaptive Wiener¯ltering, 6 curvelet domain¯ltering, 19,20 contourlet domain¯ltering, 21 Csiszars I-divergence regularization, 22 interval type II fuzzy system, 23,24 regularized image restoration based on speckle characteristics, 25 compressed sensing (CS) reconstruction, [26][27][28] sparse reconstruction, [29][30][31][32] and wavelet domain¯ltering [33][34][35][36][37] can reduce speckle noise by¯ltering in di®erent transform domain. Wavelet do-main¯ltering is widely accepted as a promising method in de-noising for OCT images.
The fundamental principle of wavelet threshold algorithm is described as follows: An estimated wavelet threshold is used to decide whether the wavelet coe±cients are noise or not. Large amplitude wavelet coe±cients are regarded as the carrier of useful information, while small amplitude wavelet coe±cients are more likely to be the carrier of noise. 36 Conventional wavelet de-noising algorithm assumes that the level of noise intensity is the same at di®erent sub-bands, so this algorithm uses a global and constant threshold for di®erent wavelet sub-bands. In fact, the speckle noise pattern is uncorrelated and di®erent in OCT images 33 and speckle noise can also be the carrier of useful information. 38 In order to protect useful information from being removed in wavelet¯ltering, it is important to evaluate suitable thresholds for di®erent wavelet sub-bands. Earlier, hard threshold and soft threshold function used in noise reduction would cause ringing e®ects 39 and Gibbs phenomenon 21,33 due to discontinuity and constant deviation in these two functions. So, it is also important to improve the wavelet threshold function.
In this paper, we propose a modi¯ed hierarchical threshold-selected algorithm and an improved wavelet threshold function based on the wavelet threshold de-noising algorithm proposed by Donoho and Johnstone. 40 The modi¯ed hierarchical threshold-selected (MHTS) algorithm can provide di®erent and suitable wavelet thresholds for di®erent wavelet decomposition levels. The improved wavelet threshold function (IWTF) has the ability to balance speckle removal and sharpness degradation. The combination of these two improvements can not only increase the contrast and signal-to-noise ratio (SNR) of OCT images, but also preserve the feature of the structural image.

Theory and Principle 2.1. Modi¯ed hierarchical thresholdselected (MHTS) algorithm
The basis theory of wavelet threshold method is described as follows: Because of the low correlation or uncorrelation between the noise and useful signals, wavelet transform makes useful signals distribute in high frequency sub-bands, meanwhile the noise signals distribute in low frequency sub-bands. After wavelet decomposition, useful signals concentrate on larger amplitude wavelet coe±cients and noise signals concentrate on smaller amplitude wavelet coe±cients. 36 For OCT images, noise intensity is di®erent at di®erent spatial or frequency positions, so the selected wavelet threshold for each decomposition level should also be di®erent. In conventional wavelet threshold de-nosing method, an unchanged global threshold is assessed and then used for each decomposition level. This kind of irrationally selected threshold could cause sharpness degradation and images' vagueness by removing some useful information during the de-noise process. Therefore, it is important to modify the thresholdselected method to achieve suitable threshold for di®erent decomposition levels.
In this paper, a hierarchical threshold is put forward based on the conventional global threshold proposed by Donoho. (As shown in Eqs. (1) and (2).) The conventional global threshold (T) proposed by Donoho is de¯ned as follows: where is the standard deviation of noisy wavelet coe±cient, HH 1 is the diagonal high frequency coe±cient at the decomposition level 1, N is the length of the signal. In practice, at high level of wavelet decomposition, wavelet coe±cients of useful signal become larger and the noisy wavelet coe±cients become smaller. If the selected global threshold is too small, it would excessively process wavelet coe±cients of useful information especially in high wavelet decomposition levels, which could cause sharpness degradation. If the selected global threshold is too large, it would retain some speckle noise in OCT images after the¯ltering process. In order to trade-o® the problem that thresholds are either too large or too small, a hierarchical threshold-selected method (as shown in Eqs. (3) and (4)) is present based on the global threshold method where j is the decomposition level, T j is the threshold at j level, HH j is diagonal high frequency coe±cient, N j is the length of HH j . Compared to the original threshold selected method, the modi¯ed method has two advantages: (1) Divided by L j , the value of which is always larger than 1. With the increasing of level j, L j becomes larger. So, the assessment thresholds T j become smaller, which meets the condition that noisy wavelet coe±cients are smaller at higher decomposition levels.
(2) An alterable parameter K is used to preserve some useful information from being unnecessarily¯ltered.

Improved wavelet threshold function (IWTF)
Wavelet threshold function is used to remove the noise mainly from low frequency wavelet coe±cients. Based on two conventional threshold functions: hard threshold function and soft threshold function, an improved wavelet threshold function is proposed to remedy the defects, such as ringing e®ect and Gibbs phenomena, 21,33,39 of these two conventional functions. The hard threshold function is shown in Eq. (5) as where " w j;k is the wavelet coe±cient after threshold processing, w j;k is the original wavelet coe±cient, j is the wavelet decomposition level, k is the spatial coordinate, T is the threshold. The discontinuities at threshold AET could cause ringing e®ect and pseudo-Gibbs phenomena, it would obscure edge details of OCT images. The blue curve in Fig. 1(a) shows the plot of the hard threshold function (assuming T is 40).
The soft threshold function is shown in Eq. (6): where sgn () is a function: if the wavelet coe±cient w j;k is larger than zero, then it will become 1, otherwise it will become À1. In soft threshold process, since there is a constant deviation between the original wavelet coe±cients and the processed wavelet coe±cients, this deviation could decrease contrast and smoothness of OCT images. The green curve in Fig. 1(a) shows the plot of the soft threshold function (assuming T is 40). In order to overcome the discontinuities and constant deviations coming from two conventional threshold functions, a new improved threshold function is presented based on a curvature function.
The de¯nition of this curvature function is as follows: The variable x of the function can be any real number on X-axis, the variable y can be any real number between (À1, 1). This function is zerosymmetric, and the degree of curvature can be adjusted by parameter . Figure 1(b) shows the plot of the modi¯ed curvature function when parameter are 0.1, 0.3, 0.5, 1, 2 and 5, respectively.
Based on the curvature function, a new threshold function is de¯ned as follows where parameter is used to realize di®erent degree of shrinkage on high frequency coe±cients. In both hard and soft threshold functions, for those wavelet coe±cients, whose absolute value are smaller than T , are totally regarded as noise and set to zero. This will inevitably remove some useful information and lead to detail de¯ciency. In the new threshold function, the wavelet¯ltering process would become more reasonable by introducing two parameters. For those with jw j;k j < T , the wavelet coe±cients after threshold process would be adjusted by choosing an appropriate parameter , rather than setting to zero in previous wavelet functions. For those with jw j;k j > T , the degree of curvature is adjusted by parameter . When is 0, the improved¯ltering process is equivalent to soft threshold function. When ! 1, the improved¯ltering process is equivalent to hard threshold function. Figure 1(c) shows the diagram of improved threshold function at di®erent degrees of curvature, parameter is 0.2, 0.5, 1.0, 2.0, 5.0, 10.0 and 20.0, respectively (assuming threshold T ¼ 40, parameter is 0.2). Figure 1(d) shows the zooming up of the area pointed by black arrow in Fig. 1(c).

Procedure of improved wavelet hierarchical threshold de-noising
As shown in Fig. 2, there are¯ve basic procedures of the wavelet-based hierarchical threshold de-noising process.

Image quality metrics parameters
Four metrics, including signal-to-noise ratio (SNR), contrast-to-noise ratio (CNR), equivalent number of looks (ENL) and XCOR, are used to assess the image quality after¯ltering process. 25 where I is the pixel value of the image on linear intensity scale, and n is the noise standard deviation of the whole image. The de¯nition of CNR 25 is shown in Eq. (13), where R is the number of ROIs choosen from the OCT images, r is the mean pixel value and r is the pixel standard deviation in the rth ROI, b and b are the mean pixel and standard deviation of the background region of the image, respectively. The CNR indicates the contrast between image feature and an area of background noise. The de¯nition of ENL 25 is shown in Eq. (14), where H is the number of ROIs, h is the mean pixel value and h is the pixel standard deviation in the h 0 th ROI. The ENL indicates smoothness in areas that should have a homogeneous appearance but are corrupted by speckle.
where I before ði; jÞ and I after ði; jÞ are the intensity value of each pixel of noisy OCT images and noisefree images, respectively. The value of XCOR is smaller than 1 and will approach 1 when the denoised images resembles the original image.

Experimental Results and Analysis
In order to demonstrate the e®ectiveness of the hierarchical threshold-selected method and improved wavelet threshold function, this technique is applied to a human¯nger skin image from a Swept-source OCT. The image size is 1024*656 pixels. The SSOCT utilizes a broadband swept source (Santec. Inc, HSL-20, 0 ¼ 1310 nm, Á ¼ 107 nm) with 12 m axial resolution and 13 m transverse resolutions. The wavelet decomposition level is 3 and processing strategy is the db6 in Daubechies group. This algorithm is realized in Matlab environment with 320.9 ms computation time for each 2D image (1024*656 pixels). The computer we used has an i7 processor. Five regions of interest (ROIs) are marked out (shown in Fig. 6(a)) in original image, that are used to calculate four metrics after the¯ltering process. Region 1 is background noise region used to estimate noisy level of the whole image. The CNR and ENL are calculated on the average of four ROIs from region 2 to 5 (marked out in the noisy image). Table 1 shows e®ect of parameter K in preserving some useful information from being excessivelȳ ltered. A 3*3 Wiener¯lter is applied to the image as a comparison. As shown in Table 1, when K ! 0:8, the CNR and ENL of MHTS are signi¯cantly larger than that of Wiener¯lter. The SNR of MHES is larger than that of Wiener¯lter when K ! 1:0 and smaller when K < 1:0. As for the XCOR, the result is another way around. Although the XCOR is decreased 0.15% to 0.35%, it is still important to get a suitable K in order to prevent some signals from being removed. From Eqs. (3) and (4), it is obvious that when K < 1, threshold T becomes smaller; when K > 1, then threshold becomes bigger. By adjusting parameter K, meanwhile and are constant at 0.2 and 5.0, respectively. As shown in Fig. 3(b), the SNR increases maximum at 9.58 dB with 3.81% sharpness degradation (when K is 1.2) compared to the original image. Meanwhile, from Fig. 3(a), CNR and ENL are also improved with the increase of parameter K. It indicates speckles are not always noisy, sometime they could also be the carrier of signals. 38 Appropriately increasing threshold T could protect these signals from being removed in wavelet¯lter.  Table 2 shows e®ect of parameter in trade-o® between speckle removing and sharpness degradation (when ¼ 5:0, K ¼ 1:2). With the decreasing of parameter from 0.5 to 0.0, the SNR has improved from 6.15 dB to 11.02 dB compared to original image, as well as showing prominent advance in both CNR and ENL (as shown in Fig. 4(a)). But it causes sharpness degradation from 2.32% to 5.19%. Figure 4(b) shows the SNR improvement and sharpness degradation for an SSOCT image, when is varied from 0.0 to 0.5. In order to balance the SNR and XCOR, a suitable parameter ¼ 0:2 is selected, the improved wavelet method can achieve  9.58 dB improvement in SNR, with sharpness degraded by 3.81%.
As it is shown in Table 3, when parameter varies from 0.1 to 5.0, SNR increases from 66.11 dB to 69.12 dB and XCOR increases from 96.11% to 96.19%. But other two metrics (both CNR and ENL) decrease. This kind of changing trend is more visible from Fig. 5, both SNR and XCOR increase at the cost of CNR and ENL decrease. It shows the ability of parameter in trade-o® between the increasing of two metrics -SNR and XCOR, and the decreasing of other two metrics -CNR and ENL.
It can be seen from the three tables that the de-noising results of the improved wavelet method is better than Wiener¯lter, soft threshold and hard threshold wavelet¯lter. By adjusting the parameters , and K, the improved threshold method has successfully balanced speckle removal and sharpness degradation, as well as protected the balance between image quality improvement and edge detail information. When ¼ 0:2, ¼ 5:0 and K ¼ 1:2, the improved wavelet method can achieve 9.58 dB improvement in SNR, with sharpness degraded by 3.81%. Figure 6 shows the comparison of an OCT¯nger image before and after di®erent¯ltering processes. Figure 6(a) is the original OCT image. Figure 6(e) is the image after improved threshold wavelet¯lter (when ¼ 0:2, ¼ 5:0 and K ¼ 1:2). The noisefree image in Fig. 6(e) shows the ability of the improved wavelet method to reduce speckle noise in OCT images. For comparison, the de-speckled images by Wiener¯ler, soft-threshold wavelet¯lter and hard-threshold wavelet¯lter are also shown in Figs. 6(b)-6(d), respectively. In order to improve the visibility of the wavelet¯lter e®ect, the 5th ROI are zoomed up (shown in Figs. 6(f) and 6(g)). Figure 6(f) is the zooming up from the 5th ROI from the original noisy image. Figure 6(g) is the same region after improved wavelet threshold¯lter. There are a lot of grains in Fig. 6(f), as a result of the existence of speckle noise. While in Fig. 6(g), the image becomes more smooth due to the speckle reducing. The improved wavelet method is not only better than those three¯lters in aspect of speckle reducing, but also capable in balance between image quality improvement and edge details preserving.

Conclusions
In OCT systems, the speckle is not only noise but also the carrier of useful information. In this paper, a modi¯ed hierarchical threshold-selected algorithm is proposed to protect signals from being removed by adjusting parameter K, which is helpful in choosing suitable thresholds for di®erent noise level. An improved wavelet threshold function is used to trade-o® between speckle removal and sharpness degradation. This improved threshold function has good continuity and adaptability, parameters and can be adjusted according to the actual feature of the OCT images. Finally, when ¼ 0:2, ¼ 5:0 and K ¼ 1:2, the improved wavelet¯lter can achieve 9.58 dB improvement in SNR, with sharpness degraded by 3.81%. The results of OCT images show that the improved threshold¯lter obtains better objective evaluation metrics and better visual discerniblity.
Except for OCT imaging, wavelet¯lters have already been applied to some other imaging systems, especially ultrasound imaging. In OCT images, the speckle noise mainly comes from lowcoherence interferometry or the imaging system itself. So, the improved wavelet we present in this paper is suitable for most OCT images/samples and its e®ectiveness is not limited to images with certain characteristics. In this paper, the algorithm is realized in Matlab. It may take more time if we transfer this algorithm to Cþþ language. But we believe with the help of graphics processing unit (GPU), we can accelerate the¯ltering process. It can be possible for real-time 2D or 3D OCT imaging.

Con°ict of Interest
The authors declare that there are no con°icts of interest related to this paper.