SAR Image Despeckling with Adaptive Multiscale Products Based on Directionlet Transform

Abstract

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 transform on the logarithmically transformed SAR images and multiply the coefficients at adjacent scales to enhance the details of image under consideration. Then, different from traditional thresholding methods, a threshold is applied to the multiscale 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.

1. 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–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–9]. Wavelet-based denoising methods usually either recover the noise-free 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–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.

2. 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 𝐼𝑥,𝑦=𝑆𝑥,𝑦⋅𝜂𝑥,𝑦. (1)

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 𝑓𝑥,𝑦=𝑔𝑥,𝑦+𝜖𝑥,𝑦, (2) 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 log-normal random variable [20] can be generated using 𝜂_log-normal=exp𝜂_normal𝑀2log𝑚+ln𝑚, (3) 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].

3. Undecimated Directionlet Transform

3.1. 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(𝑀_Λ,𝑛₁,𝑛₂)) 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 𝐝₁ (transform direction) and 𝐝₂ (alignment direction), and the coefficients are also integers. The lattice can be represented by a generator matrix 𝑀_Λ=𝑎₁𝑏₁𝑎₂𝑏₂=𝐝₁𝐝₂, (4) where 𝑎₁, 𝑏₁, 𝑎₂, and 𝑏₂∈𝑍.

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.

According to lattice theory [23], any cubic lattice ℤ² can be partitioned into |det(𝑀_Λ)| cosets of the lattice Λ, where each coset is determined by the shift vectors 𝐯_𝑘, for 𝑘=0,1,…,|det(𝑀_Λ)|−1. When the 1-D WT is applied on the pixels along the transform direction determined by 𝐝₁, after subsampling, the retained points belong to sublattice Λ^′ (Λ^′⊂Λ) corresponding to a generator matrix given by 𝑀_Λ′=2𝐝₁𝐝₂^𝑇. Such a subsampling operation allows for alignment of the retained pixels in the direction determined by the second vector 𝐝₂ 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(𝑀_Λ,𝑛₁,𝑛₂) has 𝑛₁ and 𝑛₂ 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.

3.2. 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).

4. 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.

4.1. 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] 𝑊₂^𝑗𝑓(𝑥)≤𝐾2^𝑗𝛼, (5) 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 𝑃_𝑗𝑓(𝑥)=^𝑘2_𝑖=𝑘1𝑊_𝑗+𝑖𝑓(𝑥), (6) where 𝑘₁ and 𝑘₂ 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 𝑘₁=0 and let 𝑘₂=1; then, the DWT multiscale products are 𝑃_𝑗𝑓(𝑥)=𝑊_𝑗𝑓(𝑥)⋅𝑊_𝑗+1𝑓(𝑥). (7)

Obviously, the number of the multiscale products varies with the type of transform. For UDT, the multiscale products have 2^𝑛1+𝑛2−1 components 𝑃^𝑑_𝑗𝑓𝑥,𝑦=𝑊^𝑑_𝑗𝑓𝑥,𝑦⋅𝑊^𝑑_𝑗+1𝑓𝑥,𝑦, (8) where 𝑑 represents the subband direction.

4.2. 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.(1)Apply the undecimated version of S-AWT(𝑀_𝑖,2,1), for 𝑖=1,2,…,𝑁, to the logarithmically transformed noisy SAR image 𝑓(𝑥,𝑦) up to 𝐽 scales.(2)For each coset,(i)multiply the directionlet coefficients at adjacent two scales to obtain the multiscale products 𝑃^𝑑_𝑗𝑓(𝑥,𝑦)=𝑊^𝑑_𝑗𝑓(𝑥,𝑦)⋅𝑊^𝑑_𝑗+1𝑓(𝑥,𝑦), 𝑗=1,2,…,𝐽, 𝑑=LHL,LHH,HHL,LLH,HLL,HLH,HHH (L stands for lowpass filtering and H stands for highpass filtering);(ii)compute the threshold 𝑇^𝑑_𝑗𝑝 and apply it to 𝑃^𝑑_𝑗𝑓(𝑥,𝑦) to identify 𝑊^𝑑_𝑗𝑓(𝑥,𝑦) by 𝑊^𝑑_𝑗𝑓=𝑊𝑥,𝑦_𝑑𝑗𝑓−𝑇𝑥,𝑦^𝑑𝑗𝑝𝑊^𝑑_𝑗𝑓𝑥,𝑦,𝑃^𝑑_𝑗𝑓𝑥,𝑦≥𝑇^𝑑_𝑗𝑝0,𝑃^𝑑_𝑗𝑓𝑥,𝑦<𝑇^𝑑_𝑗𝑝. (9)(3)Do the inverse S-AWT using 𝑊^𝑑_𝑗𝑓(𝑥,𝑦) to reconstruct the 𝑓_𝑖(𝑥,𝑦); then perform an exponential transform and obtain the denoised images 𝐼_𝑖(𝑥,𝑦).(4)Compute the final despeckled image 𝐼(̂𝑥,𝑦) using 𝐼(̂𝑥,𝑦)=1/𝑁∑^𝑁_𝑖=1𝐼_𝑖(𝑥,𝑦).

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.

4.3. Estimation of Threshold and Noise Standard Deviation

To derive an appropriate threshold, it is absolutely necessary to study the statistical characteristics of the multiscale products before making an estimate. Recall that the directionlet coefficients of speckle noise have been assumed to be Gaussian distribution with zero mean in the preceding section. Let 𝑊^𝑑_𝑗𝜖(𝑥,𝑦) and 𝑊^𝑑_𝑗+1𝜖(𝑥,𝑦) be zero-mean jointly Gaussian distribution with covariance matrix 𝜎𝐶=^2𝑗𝜌_𝑗,𝑗+1𝜎_𝑗𝜎_𝑗+1𝜌_𝑗,𝑗+1𝜎_𝑗𝜎_𝑗+1𝜎²_𝑗+1, (10) where 𝜌_𝑗,𝑗+1 is the correlation coefficient between 𝑊^𝑑_𝑗𝜖(𝑥,𝑦) and 𝑊^𝑑_𝑗+1𝜖(𝑥,𝑦). For the coefficients of directionlets, 𝜌_𝑗,𝑗+1 is calculated by 𝜌_𝑗,𝑗+1=∫∫𝜓_𝑗(𝑥,𝑦)⋅𝜓_𝑗+1(𝑥,𝑦)𝑑𝑥𝑑𝑦/∬(𝜓_𝑗(𝑥,𝑦))²𝑑𝑥𝑑𝑦⋅∬(𝜓_𝑗+1(𝑥,𝑦))²𝑑𝑥𝑑𝑦, where 𝜓_𝑗(𝑥,𝑦) and 𝜓_𝑗+1(𝑥,𝑦) are the directionlets used at 𝑗 and 𝑗+1 levels, respectively. The product 𝑍^𝑑_𝑗𝜖(𝑥,𝑦)=𝑊^𝑑_𝑗𝜖(𝑥,𝑦)⋅𝑊^𝑑_𝑗+1𝜖(𝑥,𝑦) has the following PDF [29]: 𝑝(𝑧)=1𝜋𝜎_𝑗𝜎_𝑗+1𝜌exp^𝑗,𝑗+1𝑥𝜎_𝑗𝜎_𝑗+11−𝜌²_𝑗,𝑗+1×𝐾₀|𝑧|𝜎_𝑗𝜎_𝑗+11−𝜌²_𝑗,𝑗+1, (11) where 𝐾₀ is the zero-th order modified Bessel function of the second kind.

The standard deviation of 𝑍^𝑑_𝑗𝜖(𝑥,𝑦) is 𝜎_𝑧=1+2𝜌²_𝑗,𝑗+1⋅𝜎_𝑗𝜎_𝑗+1. By computing the values of probability 𝑃_𝑗(𝑐)=𝑃{𝑧≤𝑐⋅𝜎_𝑧}, where constant 𝑐 varies from 1 to 5 by step length 1, the probability 𝑃_𝑗(𝑐) will be very close to 1 when 𝑐≥5 [27]. It implies that 5𝜎_𝑧 will suppress most of the data in 𝑍^𝑑_𝑗𝜖(𝑥,𝑦). In the present work, we use 𝑇^𝑑_𝑗𝑝=51+2𝜌²_𝑗,𝑗+1⋅𝜎_𝑗𝜎_𝑗+1 (12) as the multiscale products threshold to distinguish signal from noise.

As to the determination of noise level, the standard deviation of additive noise is generally estimated by 𝜎=MAD(|𝑑(𝑥,𝑦)|)/0.6745, where 𝑑(𝑥,𝑦)∈HHH₁ 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 {𝑌₁,𝑌₂,…,𝑌_𝑁} of a sequence 𝑌=(𝑌_𝑘)_𝑛∈ℕ 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) Compute 𝑛_min according to ℎ=1/4𝑁(1−𝑄) and 𝑛_min=𝑁2−ℎ𝑁, (13) and set 𝜆=√2Γ((𝑑+1)/2)/Γ(𝑑/2). (3) Denote by 𝑌₍₁₎,𝑌₍₂₎,…,𝑌_(𝑁) the sequence of observations 𝑌₁,𝑌₂,…,𝑌_𝑁 sorted by increasing norm, and compute 𝜉(𝜌)=cosh⁻¹(exp(𝜌²/2))=(1/2)𝜌+(1/𝜌)log(1+1−exp(−𝜌²)) if 𝑑=1 and 𝜉(𝜌)=(1/𝜌)𝐼₀⁻¹(exp(𝜌²/2)) if 𝑑=2, where 𝐼₀ is the zero-th order modified Bessel function of the first kind. (4) Search a smallest integer 𝑛 in {𝑛_min,…,𝑁} such that 𝑌_(𝑛)≤𝑀^∗_{{‖𝑌1‖,‖𝑌2‖,…,‖𝑌𝑁‖}}(𝑛)𝜆𝜉𝜌<𝑌_(𝑛+1), (14) where 𝑀^∗_{{‖𝑌1‖,‖𝑌2‖,…,‖𝑌𝑁‖}}(𝑛)=1𝑛^𝑛_𝑘∑₌₁𝑌_(𝑘)if𝑛≠00if𝑛=0, (15) if there exists a such integer 𝑛; then set 𝑛^∗=𝑛. Otherwise, set 𝑛^∗=𝑛_min. (5) Compute the estimate 𝜎^∗_{{𝑌1,𝑌2,…,𝑌𝑁}} of the noise standard deviation by using 𝜎^∗_{{𝑌1,𝑌2,…,𝑌𝑁}}=𝑀^∗_{‖𝑌1‖,‖𝑌2‖,…,‖𝑌𝑁‖}𝑛^∗𝜆. (16)

5. 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).

5.1. 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 (3) and 𝐿=(𝑀/𝜎)², where 𝑀 and 𝜎 are the mean and the standard deviation of the log-normal distribution, respectively; (2) then, for a specific 𝐿, we use the calculated 𝑚 and the standard normal distribution 𝜂_normal to generate the simulated speckle.

In the WHT, WMP, Cauchy, and the proposed algorithm, the orthogonal wavelet of Daubechies’ Symmlet 4 was employed and the image was decomposed into four resolution levels. Considering the computational complexity, a set of four directionlet transforms S-AWT(𝑀_𝑖,2,1)  (𝑖=1,2) has been used in the proposed method. The corresponding generator matrices are 𝑀₁=¹⁰₀₁, 𝑀₂=¹¹₋₁₁, 𝑀₃=⁰¹₁₀, and 𝑀₄=⁻¹¹₁₁. Notice that the four transforms are oriented along 0^∘, 45^∘, 90^∘, and −45^∘, respectively. Moreover, in the process of estimating noise standard deviation, two 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=10log₁₀∑^𝑁_𝑖=1𝑆²_𝑖∑^𝑁_𝑖=1𝑆_𝑖−𝑆_𝑖², (17) 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 Γ(𝑆₁,𝑆₂)=∑^𝑁_𝑖=1𝑆_1𝑖⋅𝑆_2𝑖. The correlation 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.

Table 1 

𝑆/MSE and 𝛽 values obtained by all despeckling methods applied on the 𝑆𝑦𝑛-𝑆𝐴𝑅 image.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2 

Results obtained by despeckling the Syn-SAR image. (a) Original. (b) Speckled by 9-look noise. (c) Despeckled using WHT. (d) Despeckled using WMP. (e) Despeckled using Cauchy. (f) Despeckled using the proposed method.

5.2. 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 with size 256×256 are shown in Figures 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.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 3 

(a) Noisy Suburban image. Despeckled images using (b) WHT; (c) WMP; (d) Cauchy; (e) proposed method; (f) ratio image after WHT despeckling; (g) ratio image after WMP despeckling; (h) ratio image after Cauchy despeckling; (i) ratio image after the proposed method despeckling.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 4 

(a) Noisy Mountain image. Despeckled images using (b) WHT; (c) WMP; (d) Cauchy; (e) proposed method; (f) ratio image after WHT despeckling; (g) ratio image after WMP despeckling; (h) ratio image after Cauchy despeckling; (i) ratio image after the proposed method despeckling.

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.

Table 2 

ENL and mean of ratio image for suburban and mountain images.

6. 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.

Acknowledgments

This work is supported by the National Natural Science Foundation (NNSF) of China (nos. 60872163, 61272025, and 61370110). The authors also gratefully acknowledge the helpful comments and suggestions of the editors and reviewers, which have improved the presentation.