SAR Image Despeckling with Adaptive Multiscale Products Based on Directionlet Transform

Synthetic aperture radar (SAR) images are inherently affected by multiplicative speckle noise generated by radar coherent wave. In this paper, a new despeckling algorithm based on directionlets using multiscale products is proposed. We first take an anisotropic directionlet transformon the logarithmically transformed SAR images andmultiply the coefficients at adjacent scales to enhance the details of image under consideration.Then, different from traditional thresholdingmethods, a threshold is applied to themultiscale products of the directionlet coefficients to suppress noise. Since the multiplication amplifies the significant features of signal and dilute noise, the proposed method reduces noise effectively while preserving edge structures. Finally, we compare the performance of the proposed algorithm with other despeckling methods applied to synthetic image and real SAR images. Experimental results demonstrate the effectiveness of the proposed method in SAR images despeckling.


Introduction
Over the last two decades, there is still growing interests in SAR imaging for its importance in various applications, such as search-and-rescue, high-resolution surface mapping, automatic target recognition, and mine detection.Moreover, SAR imaging systems are capable of information acquisition under all weather conditions.However, SAR images are always corrupted by a multiplicative noise called "speckle" [1] due to coherent radiation in the process of imaging.The presence of speckle noise in SAR images reduces the detection ability of targets and makes scene analysis and understanding very difficult.Thus, the removal of the speckle is a critical step in tasks such as segmentation, detection, and classification.
Many spatial-domain adaptive despeckling algorithms [2][3][4][5][6] have been proposed in the past few years.Most of these methods model the multiplicative noise and scene with certain models, design a despeckling filter or estimator based on some criterions, and finally recover the noise-free images from the observations.For example, in the Lee filter [4], the multiplicative model is first approximated by a linear combination of the local mean and the observed pixel.Then, the minimum mean-square error (MMSE) criterion is applied to determine the weighting constant used to construct the filter.Despite working well in SAR images, these traditional filters usually exhibit limitation in preserving sharp features of the images.
To address these drawbacks, a multiscale geometric analysis tool called wavelet transform has been proposed and widely used in image processing successfully [7][8][9].Waveletbased denoising methods usually either recover the noisefree image by considering wavelet coefficients as some models and constructing an estimator based on a criterion or apply hard or soft threshold to the single-scale wavelet coefficients directly.It has been shown that the wavelet thresholding algorithms can provide a better reduction of speckle noise comparing to spatial-domain filters [10].However, these thresholding methods also have the following three main drawbacks: (1) The standard two-dimensional (2D) discrete wavelet transform (DWT) is an isotropic transform, in which the filtering and subsampling operations are iterated with the same number of steps along both the horizontal and vertical directions at each scale.This means that it cannot capture edges and contours properly.(2) The downsampling in DWT is a time-variant translation and has difficulties in preserving discontinuities of image [9].(3) Applying the general threshold directly to the detail coefficients can not distinguish edges from noise effectively.Moreover, most of them do not take advantage of the interscale or intrascale dependency of the coefficients.
A large number of studies have been developed to address these problems.In [11], Zhou and Shui proposed directional windows based on contourlet transform by taking advantages of captured directional information of the images.They used local Wiener filtering in the contourlet domain and achieved better performance in removing noise.Argenti et al. despeckled SAR images successfully based on undecimated wavelet transform (or stationary wavelet transform (SWT)) [12].They used the Laplacian-Gaussian distribution to model the probability density function (PDF) of SWT coefficients and obtained the despeckled image by using the maximum a posteriori (MAP) criterion.In [7], Xie et al. proposed a Bayesian estimator by modeling the wavelet coefficients as Gaussian mixture density combining the Markov [13][14][15] random field modeling.However, all the works mentioned above only focus on the solution of one of these issues involved in wavelet-based methods.
In this paper, a novel SAR image despeckling algorithm is proposed based on two new mathematic tools called directionlet transform and multiscale products, which can reduce noise while preserving edge effectively.This is completely different from [16] which assumed the directionlet coefficients following a priori Cauchy model and employs a MAP estimator to remove noise.Our contributions consist of the following: (1) to distinguish signal from noise more efficiently, the proposed algorithm applies threshold to the multiscale products instead of the single-scale directionlet coefficients directly; and (2) a novel noise level estimator used to determine an adaptive threshold is adopted in this paper.
This paper is organized as follows.In Section 2, a brief introduction of speckle model is provided.The directionlet transform and undecimated directionlet transform are presented in Section 3. Section 4 presents the directionlet-based SAR despeckling scheme including multiscale products and the shrinkage function.Moreover, the estimation of threshold is also discussed.The performance of our proposed algorithm is evaluated and compared with the existing despeckling methods in Section 5. Conclusions are drawn in Section 6.

Speckle Model
In this work, we are interested in SAR intensity image model under the assumption that the speckle is fully developed.Let (, ) denote a SAR observation, and let (, ) and (, ) denote noise-free image and the corrupting multiplicative speckle noise, respectively.Then, the SAR image model can be expressed as  (, ) =  (, ) ⋅  (, ) . ( Since the statistical properties of speckle noise (, ) have been studied by Goodman [1], a large number of models were proposed to analyse the SAR images, such as negative exponential distribution, gamma distribution [17], log-normal distribution [10], and generalized gamma distribution [18].In [19], Arsenault and April have shown that when the image intensity is logarithmically transformed, the speckle noise is approximately Gaussian additive noise, and it tends to a normal probability much faster than the intensity distribution.Thus, it is reasonable to model the log-transformed speckle noise as a Gaussian additive noise.
Taking the logarithm of both sides of (1) yields where (, ), (, ), and (, ) signify the log values of (, ), (, ), and (, ), respectively.To ensure that (, ) follows a Gaussian distribution, we use the log-normal distribution as the speckle noise model in this work.A lognormal random variable [20] can be generated using where  and  are the mean and the median values of the distribution, respectively, and  normal is a standard normal random variable.The equivalence between  (number of looks) in a speckle image and  in (3) has been given in [10].

Directionlet Transform.
As it is well known, the anisotropic wavelet transform (AWT) whose frequency decomposition is shown in Figure 1 outperforms the standard 2D DWT in the representation of edges and contours.However, the filtering and subsampling of AWT are only along the horizontal and vertical directions, which limits the capability of representation of oriented contours to some extent.To improve the performance of representation, a skewed anisotropic wavelet transform (S-AWT( Λ ,  1 ,  2 )) called directionlet transform, which is based on integer lattices, was proposed by Velisavljevic [21,22].A full-rank integer lattice Λ consists of the points obtained as linear combinations of two linearly independent vectors denoted as d 1 (transform direction) and d 2 (alignment direction), and the coefficients are also integers.The lattice can be represented by a generator matrix where  1 ,  1 ,  2 , and  2 ∈ .
According to lattice theory [23], any cubic lattice Z 2 can be partitioned into | det( Λ )| cosets of the lattice Λ, where each coset is determined by the shift vectors v  , for  = 0, 1, . . ., | det( Λ )| − 1.When the 1-D WT is applied on the pixels along the transform direction determined by d 1 , after subsampling, the retained points belong to sublattice Λ  (Λ  ⊂ Λ) corresponding to a generator matrix given by .Such a subsampling operation allows for alignment of the retained pixels in the direction determined by the second vector d 2 and efficient iteration of transform steps.Notice that if the lattice Λ is partitioned into more than one coset, the filtering and subsampling should be performed in each coset separately.The directionlet transform obtained as a combination of the lattice-based filtering and subsampling and the frequency decomposition used in the AWT results in the skewed AWT (S-AWT).Recall that, the S-AWT( Λ ,  1 ,  2 ) has  1 and  2 transforms in one iteration step along the transform and alignment directions, respectively, and the skewed transforms are applied to all cosets of the lattice Λ separately.The basic functions of the S-AWT are called directionlets since they are anisotropic and have a specific direction.

Undecimated Directionlet Transform. It has been shown
that undecimated transforms led to a better performance than critically sampled transforms [24], especially in the application of signal denoising.For this reason, undecimated directionlets transform (UDT) is employed in our proposed despeckling algorithm.UDT is also an anisotropic transform based on lattice, but it is designed to overcome the lack of translation-invariance of the directionlet transform.It is carried out by removing the downsamplers and upsamplers of the 1-D WT in directionlet transform and upsampling the filter coefficients by a factor of 2 −1 at the th level.In other words, a zero is placed between the adjacent coefficients of   (  ) at th level to form the next level filter  +1 ( +1 ).

The Proposed Despeckling Algorithm Using Multiscale Products
In this section, our goal is the design of dirctionlet-based despeckling algorithm for SAR images using multiscale products.To design the despeckling algorithm, we first discuss the multiscale products of directionlet transform.

Multiscale Products.
Signal and noise have totally different behaviors in wavelet domain.This behavior can be analyzed by using the mathematical concept of Lipschitz regularity [25].The relation between wavelet coefficients and Lipschitz component satisfies [26] where  is the decomposition scale,  is a positive constant, and  is the Lipschitz component.Equation ( 5) implies that if the uniform Lipschitz regularity  is positive, the amplitudes of the wavelet coefficients should increase with the increasing scale.On the contrary, wavelet transform magnitudes should decrease for negative  with the increasing scale.Thus, multiplication of the DWT coefficients between adjacent scales can lead to enhancement of edge structures while weakening noise.In this paper, the multiscale products are defined as where  1 and  2 are nonnegative integers.
In [27], Bao and Zhang have pointed out that it is sufficient to implement the multiplication of two adjacent scales in practice; thus, we let  1 = 0 and let  2 = 1; then, the DWT multiscale products are Obviously, the number of the multiscale products varies with the type of transform.For UDT, the multiscale products have 2  1 + 2 − 1 components where  represents the subband direction.

Despeckling Algorithm.
In this subsection, we present the complete despeckling algorithm based on multiscale products.It has been known that denoising with hard threshold sometimes exhibits oscillation (i.e., pseudo-Gibbs phenomenon) in the vicinity of discontinuities.Although soft-threshold denoising produces less oscillation than hard threshold, it is more prone to blurring the edges of image.Thus, to get satisfied results, we employ a garrote shrinkage function [28] that provides a good compromise between the hard and the soft shrinkage functions to eliminate noise.The proposed algorithm is summarized as follows.
It can be seen that the performance of despeckling algorithm depends crucially on the threshold value    .Thus, the main task of the rest is how to estimate     properly.
As to the determination of noise level, the standard deviation of additive noise is generally estimated by  = MAD(|(, )|)/0.6745,where (, ) ∈ HHH 1 and MAD denotes the median absolute deviation operation.Despite the good characteristics (high breakdown point, good efficiency, etc.) of the MAD, its performance degrades significantly when the proportion of outliers increases.In other words, the MAD estimator is inaccurate for those images containing massive fine structures.To address this problem, we introduce a new robust scale estimator of the noise standard deviation.This estimator is called the -dimensional adaptive trimmed estimator (DATE) [30] which automatically adapts to the observations where the signals have unknown probability distributions and unknown probabilities of presence less than or equal to one-half in presence of independent additive white Gaussian noise.Here, we briefly review and summarize the DATE algorithm as follows.(1) Input a finite subsequence sequence { 1 ,  2 , . . .,   } of a sequence  = (  ) ∈N of independent -dimensional real random vectors, a lower bound  ∈ [0, ∞) for the minimum signal-to-noise ratio (SNR), and a probability value  ≤ 1 − /( − 2) 2 .(2) Compute  min according to ℎ = 1/√4(1 − ) and and set  = √ 2Γ(( + 1)/2)/Γ(/2).(3) Denote by  (1) ,  (2) , . . .,  () the sequence of observations  1 ,  2 , . . .,   sorted by increasing norm, and compute () = , where  0 is the zero-th order modified Bessel function of the first kind.(4) Search a smallest integer  in { min , . . ., } such that where if there exists a such integer ; then set  * = .Otherwise, set  * =  min .( 5) Compute the estimate  * { 1 , 2 ,...,  } of the noise standard deviation by using

Experimental Results
In this section, experiments are taken on both synthetic images and real SAR images.We compare the performance of our proposed algorithm with other methods including hard thresholding method based on DWT (referred to as WHT), multiscale products method based on DWT (referred to as WMP), and the method used in [16] (referred to as Cauchy).

Experiment on Synthetic Data.
To study the performance of the proposed algorithm in smoothing and edge preservation, we expected the experimental speckle-free image containing different types with various contents.Thus, an aerial image (called "Syn-SAR") which was obtained by cropping "westconcordorthophoto" found in MATLAB Toolbox has been chosen.We obtained the speckled image by multiplying the originals with speckle noise modeled in (3).In our experiments, the test image was of size 256 × 256, and we considered three different levels of simulated speckle noise corresponding to  = 3,5, and 9. To simulate the speckle, we adopt the following two-step approach: (1) choose  = 1 and generate a standard normal distribution; then determine a one-to-one correspondences between  and  by combining parameters  and  involved in DATE are set as 4 and 0.95, respectively.
To evaluate the performances of these despeckling approaches, we computed the signal-to-mean squared error (/MSE) ratio to measure the quality of noise suppression.The /MSE is defined as [10]  MSE dB = 10 log 10 ( where  is the image size and  and Ŝ are the original and the denoised images, respectively.This measure corresponds to the classical SNR in the case of additive noise. In addition to the above evaluation measure, we also used a parameter  to evaluate the performance of edge preservation.It is originally defined as [31]  = Γ (Δ − Δ, Δ − Δ) √Γ (Δ − Δ, Δ − Δ) ⋅ Γ ( Δ − Δ, Δ − Δ) , (18) where Δ and Δ are the highpass-filtered versions of the original image  and the denoised image Ŝ, respectively, obtained with a window size of 3 × 3 standard approximation of the Laplacian operator.The overline operator represents the mean value, and measure, , should be close to unity for an optimal effect of edge preservation.Figure 2(a) illustrates the original test image Syn-SAR.The values of /MSE and  obtained by applying all the methods to the test image are listed in Table 1.It can be seen from the table that the proposed method provides larger values of /MSE in comparison to other methods; thus, it indicates a better ability to suppress the speckle noise.Table 1 also shows that the values of  obtained by our approach are larger than those obtained by the other three methods.It should be noted that the performance of WMP is increasingly superior to Cauchy with the increase of noise intensity; this is probably because multiscale products play a more important role in diluting noise.For a visual comparison, the speckled images ( = 9) and the corresponding despeckled images are shown in Figures 2(b)-2(f), respectively.They also indicate that the proposed method achieves better visual quality than the WHT, WMP, and Cauchy methods.Obviously, the result of visual comparison is consistent with the values of /MSE and  measures.

Experiments on Real SAR Images.
In order to further study the advantages of the proposed method, we also perform experiments on real SAR images.The wavelet used in the mentioned approaches and the initial values of  and  involved in DATE are set same as to the preceding simulations.The two noisy images Suburban and Mountain  3(a) and 4(a), respectively.Since the noise-free images are not available, the values of the /MSE and  can not be appropriately used to evaluate the improvement of our method.Fortunately, we can use the equivalent number of looks (ENL) and the ratio image (/ Î) [32] to assess the despeckling performances.It is important to note that the ENL must be calculated in a homogeneous region which has been highlighted with a rectangle in the original images and that the ratio should have the characteristics of pure speckle if the speckle is fully developed in areas.
To assess the quality of despeckled images, we show the resultant images in Figures 3(b)-3(i) and 4(b)-4(i) and list the values of ENL and mean of ratio image in Table 2. From the figures and the table, we can see that the proposed method achieves the best performance.The results seem to be consistent with the previous simulation results.We attribute the better performance of the proposed method to the application of directionlet-based multiscale products with the ability of amplifying the significant features and diluting noise and to the superior performance of garrote shrinkage function.

Conclusion
In this paper, a new SAR image despeckling method using multiscale products based on directionlets was proposed.We multiplied the adjacent scale subbands to amplify significant features and then applied the threshold to multiscale products instead of to the single-scale directionlet coefficients directly.The experimental results on test SAR image showed that the proposed method reduces speckle effectively while preserving edge structures.Due to application of directionlet transform, multiscale products, and a more complicated method for threshold estimation, the proposed despeckling algorithm requires a relatively high computational demand, but it is not a key issue in application for now.For future work, we plan to focus on how to determine a more appropriate threshold, which is important for a despeckling algorithm.

Figure 1 :
Figure 1: Frequency decomposition of AWT, where two decomposition levels are used.1D wavelet transform steps along two directions (horizontal and vertical) are not equal to each scale.

Table 1 :
/MSE and  values obtained by all despeckling methods applied on the - image.

Table 2 :
ENL and mean of ratio image for suburban and mountain images.
with size 256 × 256 are shown in Figures