1. Introduction

Optical Coherence Tomography (OCT) is a relatively novel biomedical imaging technology that allows in vivo, non-invasive, high-resolution imaging of biological tissue at superficial depths of 1–2 mm [1]. OCT is based on low-coherence interferometry, which utilizes the spatial and temporal coherence properties of optical waves backscattered from biological tissue [2]. As any imaging technique that is based on detection of coherent waves, OCT images are subjected to presence of speckle. Speckle in OCT tomograms is dependent both on the wavelength of the imaging beam and the structural details of the imaged object [3]. Therefore, speckle carries information both about the morphology of the imaged object as well as a noise component, and the latter is responsible for the granular appearance of OCT images. Speckle noise can obscure small or low reflectivity features in OCT images, thus degrading the image quality of OCT tomograms. Furthermore, it can impede or limit the performance of image segmentation and pattern recognition algorithms that are used to extract, analyze, and recognize diagnostically relevant features.

Development of successful speckle noise reduction algorithms for OCT is particularly challenging, because of the difficulty in separating the noise and the information components of a speckle pattern. In the past, extensive research has been conducted both in the fields of medical imaging and remote sensing for suppressing speckle noise. A number of digital speckle reduction algorithms such as adaptive filters [3] and rotating kernel transformation [4] have been developed or adapted specifically for improving image quality of OCT tomograms. Filtering techniques based on the rotating kernel transform can produce good contrast enhancement of image features, however they result in significant edge blurring when strong noise reduction is required [5]. A wavelet based soft thresholding technique is described in [5]. It computes the undecimated wavelet transform and applies soft thresholding to the horizontal, vertical, and diagonal subbands. The threshold is obtained using the statistics of the wavelet coefficients. The wavelet based technique described in [5] does not reduce the image sharpness significantly but the execution time for the algorithm is about 7 min using Matlab implementation. In a different study [6], the performance of various digital filters was compared by applying the algorithms to an OCT tomogram of a bovine retina acquired ex vivo. More specifically, the results presented in reference [6] indicate that non-orthogonal wavelet-transform-based filters together with enhanced Lee and adaptive Wiener filters can reduce speckle and increase SNR while preserving fine details of image. The nonorthogonal wavelet-transformed filter applies a separate threshold for each individual wavelet level. Enhanced Lee filter is an improved version of the original Lee filter which a minimization of the mean square error estimate of the true image with the multiplicative noise model transformed into an independent additive noise [7]. The Adaptive Wiener filter is also determined from the minimization of the mean square error estimate of the true image. It is implemented based on the local statistics (mean and variance) within the smoothing window centered on each pixel of the input image [6].

Most of the algorithms are limited in the amount of speckle that can be reduced because of their complexity. In addition to algorithms, modifications to the OCT system design have been suggested to reduce speckle noise. One of the techniques uses a multimode fiber and a spatially coherent broadband light source to reduce speckle noise [8]. Spatial and frequency compounding are alternative techniques for reducing speckle noise in OCT images, however, these techniques require significant modifications to design of the OCT system. Recently, a spatial diversity technique based on shifting the focal plane of the probe beam was introduced to suppress speckle noise [9]. In a different study, angular compounding technique based on path length encoding was presented to reduce presence of speckle in OCT images [10]. Speckle noise reduction in OCT images was recently achieved by use of angular compounding at multiple backscattering angles [11]. The study described in reference [6] proves that the combination of experimental techniques (specifically angular compounding) and digital filter processing can result in overall superior image quality enhancement of OCT tomograms.

In this paper, we present a robust, spatially adaptive speckle noise reduction technique in the wavelet domain, based on the interval type II fuzzy sets (also called ultra-fuzzy sets). This method utilizes soft thresholding to suppress the wavelet coefficients. The approach is a modification and improvement to previously published work [12] that uses type I fuzzy sets for additive Gaussian noise reduction.

2. Theory

In 1989 Stephane G. Mallat, introduced a multi-resolution representation of an image that can be used to analysis the information content present in both coarse and fine details. Wavelet transform is a multi-resolution decomposition tool, which allows an image to be represented at various resolutions or levels. By filtering and scaling the input image, the image space is divided into sequences of four lower resolution (or lower scale) components containing the approximation, horizontal, vertical, and diagonal detail coefficients represented by W_φ(j,m,n), W^H_ψ(j,m,n), W^V_ψ(j,m,n), and W^Dψ(j,m,n) respectively. The variable j represents the decomposition scale level and(m,n) represents the image spatial location. The crucial image information is compressed and represented along several resolution scales by few of these large valued coefficients [13]. Thus, by analyzing images at various resolutions, one can remove the unwanted noise, represented by the low valued coefficients.

In image processing, the application of a wavelet based de-noising method consists of the following three steps: i) computing the two-dimensional discrete wavelet transform (2D-DWT) of the image at various levels of decomposition ii) removing noise from the detail wavelet coefficients (W^H_ψ, W^V_ψ, and W^Dψ) by soft thresholding iii) finally, reconstructing the enhanced image using inverse 2D-DWT. For an image, Mallat’s 2D-DWT can be implemented using 4 frequency channel filter banks corresponding to the four orientations. The wavelet transform of natural images has special properties like interscale and intrascale dependencies and as a result de-noising using wavelet transform has proven to be very effective. Interscale dependency refers to large wavelet coefficients propagating across the scales and intrascale dependency refers to neighbourhood of similar valued coefficients (large or small) at each scale.

In OCT images, speckle noise can be approximated as a multiplicative noise such that,

where s(m,n) represents the noise free OCT image, f(m,n) is the noise observation of s(m,n), n(m,n) and n_a(m,n) are multiplicative speckle and additive noise respectively, and(m,n) is the variable of spatial locations. In Eq. (1), the additive noise component comes from the shot noise, light intensity noise, and electronic noise inherent in OCT imaging system and it can be ignored because it is significantly small compared to the multiplicative speckle noise [14]. For developing speckle reduction filters in the wavelet domain, multiplicative noise has to be transformed into additive noise. This can be achieved by use of logarithmic transformation, resulting in Eq. (2).

Once, the image has been transformed into logarithmic scale, the 2D-DWT can be applied:

In Eq. (3), f_j,d denotes the noisy wavelet coefficient of the observed image at scale j and at orientation d, s_j,d is the noise-free wavelet coefficient of the image and n_j,d(m,n)represents the wavelet coefficient of the speckle. Additive noise in the image domain remains additive in the transformed wavelet domain due to the linearity of the wavelet transform. The next step in the wavelet based de-noising algorithm is thresholding the detail wavelet coefficients, f_j,d(m,n)where d ∊{H,V,D}, by use of a soft threshold. In soft thresholding if the magnitudes of the coefficients are below the threshold then they are set to zero while the coefficients with the magnitudes above the threshold are scaled towards zero. This is because the coefficients that contain mostly noise should be reduced to negligible values, while the ones containing a substantial noise free component should be reduced less. In wavelet domain different types of noise are associated with small magnitude coefficients. Important image structures are contained within the magnitude of the high coefficients. The coefficients around the threshold contain both noise and image features of interest. Therefore, an optimal threshold is reached when most of the coefficients bellow it are noise and the coefficients above it represents image features of interest.

Due to the ambiguity of choosing a suitable threshold, it is appropriate to apply the fuzzy set theory in such situations. Fuzzy set theory is a soft computing technique that deals with the imprecision and vagueness of human understanding systems. It is able to directly model and minimize the effect of uncertainty in the calculated threshold. Fuzzy set theory was first introduced by L. A. Zadeh in 1965 [15] and since then has been used extensively in image processing. Using fuzzy set theory, more specifically, using fuzzy rule-based approach, the wavelet thresholding can be expressed as a fuzzy wavelet thresholding. The fuzzy wavelet thresholding procedure can be divided into five main steps which are discussed in detail below.

Step 1. Establish the image features, each corresponding to linguistic labels. Each feature is then fuzzified by their appropriate membership function (MF). Linguistic label is a verbal scale of qualitative linguistic estimates [16], for expressing one’s own subjective judgments about both the support and the negation of a proposition [16]. For example, “likely in very high degree” and “unlikely in high degree” are two examples of linguistic labels that can be used to describe the proposition “Is that person tall?” These linguistic labels help to translate human reasoning to digital operations. These labels are represented by fuzzy sets in fuzzy set theory. The MF is a mathematical function that specifies the degree to which a given input belongs to the set defined by the linguistic labels. The output of the MF is between 0 and 1. In this manuscript, two fuzzy variables are defined and the following describes it.

At position (m,n), if the wavelet coefficient has a magnitude that is large, then the coefficient at (m,n) represents image feature of interest for “almost certain” and should not be set to zero but scaled towards zero. The reason for this statement is due to the property of image edges appearing as large magnitude wavelet coefficients due to the rapid changes in intensity values. Large image structures can be found in many scales of the wavelet decomposition, but small image details can only be revealed in several fine scales and to determine small image details, interscale correlation involving two adjacent scales can be used [17]. A hierarchical correlation map that considers both the correlation between neighbouring coefficients (intra-scale) and correlation between two adjacent levels (inter-scale) can be obtained to determine fine image features. Figure 1 shows the intra-scale and inter-scale relationship among the wavelet coefficients. In the correlation map that is generated, large values indicate position of edges in the original image, and zero coefficients correspond to smooth areas which are noise regions [17]. This leads to the following statement: at position (m,n), if the correlation map value is large, then the coefficient at (m,n)represents fine image structures of interest for “almost certain” and should not be set to zero but scaled towards zero. Conversely, a correlation value close to zero indicates an area which needs to be smoothed due to noise.

 

Fig. 1. Two-dimensional wavelet decomposition at level j and j+1. The intra-scale and interscale relationship among the horizontal (H), vertical (V) and diagonal (D) detail coefficients is shown.

Download Full Size | PDF

Using the linguistic labels, large magnitude wavelet coefficient at (m,n) and large correlation map value at position (m,n), two fuzzy variables can be defined as follows: Large magnitude wavelet coefficient:

Large correlation map value modified from reference [17] to account for the undecimated 2D wavelet transform:

Fine image structures which do not appear as local maxima will now be revealed in the correlation map [17]. Due to oriented edges being visible in certain subbands, subands of different orientations should be treated separately. As a result, the correlation map is found for the horizontal, vertical, and diagonal subbands separately, i.e. d∊{H,V,D}. The correlation map value x_j,d(m,n) at decomposition level j and orientation d is based on nine wavelet coefficients in a 3×3 window (N=1) at level j and orientation d, and the nine coefficients in the next level j+1 and orientation d. The subbands in each of the wavelet levels have the same size due to the Stationary Wavelet Transform (SWT) used in this manuscript. As a result the coefficients have one-to-one correspondence in terms of the position(m,n).

Now, two corresponding Interval Type II Fuzzy Sets (IT2FS) can be assigned to these fuzzy variables. Type II fuzzy set A corresponds to the large magnitude wavelet coefficient and type II fuzzy set B corresponds to the large correlation map value. The corresponding type II sigmoid MFs are presented in Fig. 2(a) and Fig. 2(b). The shaded regions in the figures correspond to the areas where the uncertainty lies in determining if the coefficient is noisy or not. The general equation for a sigmoid MF is given in Eq. (6) and it has two parameters associated with it, its center (c) and its width (w).

The upper and lower MFs are determined based on Eq. (7.1) and Eq. (7.2).

According to Tizhoosh [18], β should be between 1 and 2 and β≫2 is usually not meaningful for image data. In this manuscript, β was set to a value of 1.2. The soft threshold obtained by applying the technique in Ref. [19] was used as the center c of µA(x). This is because fuzzy set A corresponds to large wavelet coefficients and the soft thresholding procedure described in Ref. [19] utilizes this information to suppress the noise. The threshold is obtained by computing the ratio of the noise variance and the standard deviation of each of the subbands. In this manuscript, the noise variance was estimated using the background region in the vertical orientation of the first subband, in a similar way as described in Ref. [5]. Thus, it is evident, that at each scale and at each orientation, the MF of fuzzy µ_A(x) changes. The width w of µ_A(x) and µ_B(x) was set to a value of 0.3 after optimizing for SNR improvement while using the same center for µ_A(x) and µ_B(x). Finally, the location of the MF of µ_B(x) is controlled by a parameter η, which is the center, c of µ_B(x). By adjusting η the amount of speckle removed can be varied. Thus, allowing for more or less smoothing depending on the input image being used and the amount of speckle present.

 

Fig. 2. Interval Type II Fuzzy Membership functions for the fuzzy variables a) “large magnitude wavelet coefficient” and b) “large correlation map value”.

Download Full Size | PDF

In Fig. 2, the fuzzy membership value is near one for both large magnitude wavelet coefficient and large correlation map value. Similarly, if the magnitude of the wavelet coefficient or correlation map value is small or near zero, then the fuzzy membership value is near zero. However, if the magnitude of the coefficient and the correlation map value are neither large nor small, then the fuzzy membership value takes on a range of values. An example of this is depicted for both fuzzy sets in Fig. 2. This range is defined by the upper and lower membership values and it incorporates the uncertainty involved in deciding if the magnitude of the wavelet coefficient and the correlation map value near the threshold is indicating if noise is present at that location. The membership values represent the degree by which a wavelet coefficient represents an image feature of interest. High membership value corresponds to high degree that the wavelet coefficient belongs to an image feature of interest. Likewise, a low membership value corresponds to a low degree that the wavelet coefficient belongs to an image feature of interest.

Step 2. The next step in fuzzy wavelet thresholding is to set up a knowledge-base that is composed of fuzzy if-then rules. Utilizing the fuzzy variables representing the image features the following fuzzy rules are obtained:

Rule 1: If the magnitude of the wavelet coefficient at (m,n) is “large” AND the correlation map value at (m,n) is “large” then “scale towards zero the wavelet coefficient at(m,n).”

By adding new fuzzy rules, the algorithm can be changed effortlessly to include other information from the image to further improve speckle reduction performance.

Step 3. For each rule, the fuzzified inputs are combined and a rule strength is obtained. This rule strength is used as the consequence of the rule. The strength of rule 1 is computed via the algebraic product triangular norm and is given by the following:

Since, IT2FS is being used; there are two consequences for the rules each corresponding to the upper and lower limit. This consequence represents the degree by which the wavelet coefficient should be scaled.

Step 4. Once the consequences from various rules have been obtained, they are combined to obtain an output distribution range. In this manuscript only one fuzzy rule is being used and thus, this step can be ignored. In general, the probabilistic sum triangular conorm is used to combine various rules.

Step 5. Finally, a value representing the scaling factor for the wavelet coefficient at (m,n) is obtained by type reducing and defuzzifying the output from step 4. There is more than one approach to defuzzifying the output distribution; here the average of the upper and lower value is utilized.

Now, a new coefficient can be obtained at position (m,n) which has been reduced of noise. This is obtained by applying Eq. (10). The filtering is performed at each scale of the wavelet domain and for each of the three detail coefficients (d∊{H,V,D}).

Once, all the coefficients have been found, applying the inverse 2D-DWT to an input OCT tomogram will result in a de-noised final image.

The fuzzy type II speckle noise removal filter was implemented in MATLAB (v. 7.0.1) and its performance was evaluated using well-known speckle-reduction performance metrics for images; they include SNR, contrast-to-noise ratio (CNR), equivalent number of looks (ENL) and edge preservation (∊) [5, 6, and 20]. The smoothness of a homogenous region of interest is measured by ENL and CNR is a measure between an area of image feature and an area of background noise. Edge preservation is a correlation measure that indicates how the edges in the image are degrading. A value close to 1 indicates the filtered image is similar to the reference image. The image quality metrics definitions as described in references [5, 6, and 20] are presented below:

where I and σ²_n in SNR represents the linear magnitude image and the variance of the background noise region in the linear magnitude image respectively. µ_h and σ²_h in ENL represents the mean and variance of the h^th homogenous region of interests respectively. µ_b and σ²_b in CNR represents the mean and variance of the same background noise region as in SNR and µ_r and σ²_r represents the mean and variance of the r^th region of interest which includes the homogeneous regions as well. In the edge preservation measure, ΔI and ΔÎ represent the Laplace operator performed on the original image I and the filtered image Î respectively. Also, ΔĪ and ΔÎ represent the mean value in the r^th region of interest of ΔI and ΔÎ respectively.

Following a previous work [5], all the calculations were performed in logarithmic domain of an OCT image except for the SNR calculation. Six regions of interest were used for the algorithm test, including 3 which were homogenous regions (R=6,H=3). The orthogonal Daubechies wavelet, more specifically “db4”, was utilized for the 2D-DWT with 4 levels of decomposition. The reason for using the fourth-order Daubechies wavelet is it brings smoothing effects due to the longer filter lags [21] and this reduces the pixilation effect seen when lower-order Daubechies wavelet is utilized. After various experiments with different levels of decomposition, 4 levels were chosen as the optimal level. The 2D-DWT was implemented using an undecimated filter bank and this formulation is known as the SWT [22]. Although SWT requires excessive computer memory and longer CPU run-time, it is able to denoise and detect edges better [22]. If Mallat’s 2D-DWT is used then the algorithm will run at half the speed but because Mallat’s 2D-DWT implementation uses the decimated filter bank, it is prone to Gibbs artefacts caused by sharp edges in the image [5]. The digital implementation of SWT algorithm was done using the Fast Wavelet Transform (FWT) technique outlined in [23]. The image is convolved with the “db4” decomposition low-pass and high-pass filters in various sequences. The next scale of the image in the wavelet transform is obtained by up sampling the low-pass and high-pass filters by a factor of 2. As mentioned earlier in the manuscript, the four frequency channel filter banks for the 2D-DWT are the four sequences of convolution involving decomposition “db4” low-pass and high-pass filters (low-low, low-high, high-high, and high-low). The Matlab conv2() function was utilized for the convolution. Similarly, the image is reconstructed by up sampling the low scale images and convolving it with the “db4” low-pass and high-pass reconstruction filters.

The performance of the novel fuzzy type II algorithm was compared with the results obtained with an adaptive Wiener filter [24] that was provided by Matlab via the wiener2() function. Also, a modified adaptive Lee filter implemented using an algorithm described in reference [7] was applied to the image. As mentioned before, the adaptive Wiener filter and the Lee filter were previously applied to OCT images and they showed valuable results [6].

3. Results and discussion

To evaluate the performance of the new fuzzy type II algorithm, it was applied to images acquired in-vivo from a human finger tip with a high speed, high resolution Fourier Domain OCT system (FD-OCT) operating in the 1060nm wavelength region. The FD-OCT system was powered with a superluminescent diode (Superlum, λc=1020nm, Δλ=108nm, P_out=10mW), which provided 3. 5µm (axial) and 10µm (lateral) resolution in biological tissue. The FD-OCT system was equipped with a fast InGaAs linear array, 1024 pixel CCD camera (SUI, Goodrich) with a readout rate of 47kHz. The spectrometer was designed to fit the entire spectrum of the SLD, thus providing a maximum scanning range of 1.4mm in air. Two dimensional images (1000 A-scans×512 pixels) of a human finger tip were acquired with 1.3mW incident power and 97dB SNR at a rate of 46 frames/s.

The adaptive Wiener, the adaptive Lee and the fuzzy type II algorithms were applied to an OCT image of a fingertip with dimensions 1000×512(A-scans × pixels). Fig. 3A shows the unprocessed OCT fingertip tomogram that has grainy appearance due to the presence of speckle. Six regions of interest (homogeneous regions labelled 2, 4, and 6) and background region (labelled 1) above the tissue were selected and marked on the image. A region in the image, marked with a blue line box and containing what we believe are sweat glands, was enlarged by 2× and shown as an inset in Fg. 3(A). Figure 3(B) shows the same image after adaptive Wiener filtering and a similar inset shows the reduced appearance of speckle pattern. Red arrows in the image mark the locations of sweat glands in the epidermal layer of the fingertip. Figure 3(C) shows the OCT image after adaptive Lee filtering. Both the Wiener and Lee filters are based on the optimal minimum mean square error estimator of the true image. Both filters appear to reduce the speckle pattern, yet in homogenous regions some speckle patterns are still clearly visible. The type II fuzzy filtered image, shown on Fig. 3(D) exhibits clear edges of tissue features (for example in the sweat glands) and significantly reduced appearance of speckle patterns. Note that the same gray scale was used in all 3 images and the scale was not altered after the application of the Lee, Wiener and Fuzzy type II filters.

 

Fig. 3. Wavelet denoising of a human finger tip image acquired with the FD-OCT 1060nm system: A) Original OCT image. B) Adaptive Wiener filtered image. C) Adaptive Lee filtered image and D) type II fuzzy wavelet filtered image. The images are 1000×512 (lateral × axial) pixels, corresponding to 1 mm×0.8 mm physical dimensions. The red and yellow line boxes in A) mark the selected regions for the image metrics evaluation. The red arrows in B)-D) point at sweat glands in the epidermis. The blue line box in A) marks a region containing sweat glands, that has been enlarged and shown as a 2× inset in image A). Similar insets were generated for the processed images B) - D). Enlarged copies of all four insets are shown in E) - H) for close comparison of the performance of the three wavelet algorithms. E) Original unprocessed image, F) Adaptive Wiener, G) Adaptive Lee and H) Fuzzy Type II set algorithms.

Download Full Size | PDF

Enlarged views of the 4 insets are presented in Fig. 3(E) (unprocessed image), Fig. 3(F) (Adaptive Wiener filter), Fig. 3(G) (Adaptive Lee filter) and Fig. 3(H) (Fuzzy type II filter) for easy comparison of the effect of filtering on the presence of a speckle pattern. It is clear from Fig. 3(E)–3(H), that both the Adaptive Wiener and Adaptive Lee filters reduce significantly speckle noise in the images, while preserving the edges of the sweat gland. Some residual speckle pattern is still visible in the background tissue on Figs. 3(F), 3(G). The Fuzzy type II filter shows significantly improved reduction in the speckle patter while still preserving clear edges of the sweat gland.

Quantitative comparison of the filtering effect for all three algorithms was performed by evaluating the image quality metric values for all four images [Figs. 3(A)–3(D)]. The results from the quantitative analysis are presented in Table 1. The use of the Wiener and Lee filters results in some improvement in the image SNR (up to 3.6 dB and 5 dB respectively) and CNR. Image regions labelled 2, 4, and 6 in Fig. 3(A) were used for the ENL calculation and regions labelled 2 to 7 were used for the CNR calculation. The application of the Fuzzy type II filter results in more significant image quality improvement (up to 10 dB increase in SNR), which is also apparent by the change in image contrast in Fig. 2 (bottom). Clearly, the proposed Fuzzy type II based algorithm enhances the appearance of minute morphological features, in this case, the presence of sweat glands in the epidermis. Also, table 1 shows significant improvements in CNR and ENL in the wavelet processed images as compared to the original OCT image. Also, the fuzzy wavelet filtered image is close to the original image as its edge preservation value is greater than the other filtered image edge preservation value. The CPU time required to run the algorithms was recorded (table 1) and shows that the fuzzy type II is the slowest of the three algorithms, however it is much faster as compared to the algorithm suggested in Ref. [5].

Table 1. Image quality metrics for human finger tip OCT image

View Table

Optimization of the image SNR with the fuzzy type II algorithm can be achieved by varying the parameter η to vary the amount of speckle removed by the fuzzy wavelet thresholding filter. Maximum SNR of about 10dB is achieved for η=0.4325. Figure 4 shows the image SNR and a function of the parameter η.

Our technique shows that speckle noise can be significantly reduced using an interval type II fuzzy set based thresholding in the wavelet domain. The results from the experiment verify that interval type II fuzzy sets can reduce speckle pattern more effectively compared to standard filters. The Matlab implementation of the algorithm took approximately 40 seconds to execute on a computer with Pentium IV 3.0 GHz processor with 1.0 GB of RAM for test image dimensions of 1000×512 (lateral × axial) pixels. The algorithm written in Matlab was optimized to use matrix manipulation by avoiding the use of nested loops. The wavelet transform was implemented with the FWT method from [23] with the help of Matlab built in convolution functions. The proposed speckle reduction algorithm requires significantly less run time as compared to previously published algorithms, while resulting in greater improvement in the image SNR [5]. The proposed algorithm requires a one-time optimization of the η parameter each time the OCT system is set up. This is a somewhat time consuming task however, once optimization is achieved the algorithm can be applied to a batch of images with a very short image processing duration.

 

Fig. 4. SNR improvement as a function of η. Maximum SNR is obtained at η=0.4325.

Download Full Size | PDF

4. Conclusion

We have developed a novel algorithm for reducing speckle noise in OCT images, based on interval Type II fuzzy thresholding in the wavelet domain. The performance of the algorithm in terms of improvement in SNR, CNR, ENL and Edginess is superior to previously published speckle reduction algorithms that have proven to provide excellent results (Adaptive Wiener and Enhanced Lee algorithms). The fuzzy type II technique relatively slower as compared to fast processing algorithms such as Adaptive Wiener, filter, however, it is much faster as compared to other image processing techniques. We demonstrated the effectiveness of the algorithm in reducing speckle noise by applying it to an OCT image of a finger tip and achieving SNR improvement of 10dB.