# SAR Image Despeckling with Adaptive Multiscale Products Based on Directionlet Transform

Authors: Yixiang Lu; Qingwei Gao; Dong Sun; Dexiang Zhang
Journal: Mathematical Problems in Engineering
 (2013)

---

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


---

## Body

## 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 Section2, 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. LetI(x,y) denote a SAR observation, and let S(x,y) and η(x,y) denote noise-free image and the corrupting multiplicative speckle noise, respectively. Then, the SAR image model can be expressed as
(1)I(x,y)=S(x,y)·η(x,y).Since the statistical properties of speckle noiseη(x,y) 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)f(x,y)=g(x,y)+ϵ(x,y),
where f(x,y), g(x,y), and ϵ(x,y) signify the log values of I(x,y), S(x,y), and η(x,y), respectively. To ensure that ϵ(x,y) 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
(3)ηlog-normal=exp(ηnormal2log(Mm)+lnm),
where M and m are the mean and the median values of the distribution, respectively, and ηnormal is a standard normal random variable. The equivalence between L (number of looks) in a speckle image and m 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 Figure1 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(MΛ,n1,n2)) 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 d1 (transform direction) and d2 (alignment direction), and the coefficients are also integers. The lattice can be represented by a generator matrix
(4)MΛ=[a1b1a2b2]=[d1d2],
where a1, b1, a2, and b2∈Z.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 ℤ2 can be partitioned into |det(MΛ)| cosets of the lattice Λ, where each coset is determined by the shift vectors vk, for k=0,1,…,|det(MΛ)|-1. When the 1-D WT is applied on the pixels along the transform direction determined by d1, after subsampling, the retained points belong to sublattice Λ′ (Λ′⊂Λ) corresponding to a generator matrix given by MΛ′=[2d1d2]T. Such a subsampling operation allows for alignment of the retained pixels in the direction determined by the second vector d2 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(MΛ,n1,n2) has n1 and n2 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 2j-1 at the jth level. In other words, a zero is placed between the adjacent coefficients of Hj(Gj) at jth level to form the next level filter Hj+1(Gj+1).


## 3.1. Directionlet Transform

As it is well known, the anisotropic wavelet transform (AWT) whose frequency decomposition is shown in Figure1 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(MΛ,n1,n2)) 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 d1 (transform direction) and d2 (alignment direction), and the coefficients are also integers. The lattice can be represented by a generator matrix
(4)MΛ=[a1b1a2b2]=[d1d2],
where a1, b1, a2, and b2∈Z.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 ℤ2 can be partitioned into |det(MΛ)| cosets of the lattice Λ, where each coset is determined by the shift vectors vk, for k=0,1,…,|det(MΛ)|-1. When the 1-D WT is applied on the pixels along the transform direction determined by d1, after subsampling, the retained points belong to sublattice Λ′ (Λ′⊂Λ) corresponding to a generator matrix given by MΛ′=[2d1d2]T. Such a subsampling operation allows for alignment of the retained pixels in the direction determined by the second vector d2 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(MΛ,n1,n2) has n1 and n2 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 2j-1 at the jth level. In other words, a zero is placed between the adjacent coefficients of Hj(Gj) at jth level to form the next level filter Hj+1(Gj+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]
(5)|W2jf(x)|≤K(2j)α,
where j is the decomposition scale, K 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
(6)Pjf(x)=∏i=k1k2‍Wj+if(x),
where k1 and k2 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 k1=0 and let k2=1; then, the DWT multiscale products are
(7)Pjf(x)=Wjf(x)·Wj+1f(x).Obviously, the number of the multiscale products varies with the type of transform. For UDT, the multiscale products have2n1+n2-1 components
(8)Pjdf(x,y)=Wjdf(x,y)·Wj+1df(x,y),
where d 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(Mi,2,1), for i=1,2,…,N, to the logarithmically transformed noisy SAR image f(x,y) up to J scales.(2)
For each coset,
(i)
multiply the directionlet coefficients at adjacent two scales to obtain the multiscale productsPjdf(x,y)=Wjdf(x,y)·Wj+1df(x,y), j=1,2,…,J, d=LHL,LHH,HHL,LLH,HLL,HLH,HHH (L stands for lowpass filtering and H stands for highpass filtering);(ii)
compute the thresholdTjdp and apply it to Pjdf(x,y) to identify W^jdf(x,y) by
(9)W^jdf(x,y)={Wjdf(x,y)-TjdpWjdf(x,y),Pjdf(x,y)≥Tjdp0,Pjdf(x,y)<Tjdp.(3)
Do the inverse S-AWT usingW^jdf(x,y) to reconstruct the f^i(x,y); then perform an exponential transform and obtain the denoised images I^i(x,y).(4)
Compute the final despeckled imageI^(x,y) using I^(x,y)=1/N∑i=1N‍I^i(x,y).It can be seen that the performance of despeckling algorithm depends crucially on the threshold valueTjdp. Thus, the main task of the rest is how to estimate Tjdp 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. LetWjdϵ(x,y) and Wj+1dϵ(x,y) be zero-mean jointly Gaussian distribution with covariance matrix
(10)C=[σj2ρj,j+1σjσj+1ρj,j+1σjσj+1σj+12],
where ρj,j+1 is the correlation coefficient between Wjdϵ(x,y) and Wj+1dϵ(x,y). For the coefficients of directionlets, ρj,j+1 is calculated by ρj,j+1=∫‍∫‍ψj(x,y)·ψj+1(x,y)dxdy/∬(ψj(x,y))2dxdy·∬(ψj+1(x,y))2dxdy, where ψj(x,y) and ψj+1(x,y) are the directionlets used at j and j+1 levels, respectively. The product Zjdϵ(x,y)=Wjdϵ(x,y)·Wj+1dϵ(x,y) has the following PDF [29]:
(11)p(z)=1πσjσj+1exp(ρj,j+1xσjσj+1(1-ρj,j+12))×K0(|z|σjσj+1(1-ρj,j+12)),
where K0 is the zero-th order modified Bessel function of the second kind.The standard deviation ofZjdϵ(x,y) is σz=1+2ρj,j+12·σjσj+1. By computing the values of probability Pj(c)=P{z≤c·σz}, where constant c varies from 1 to 5 by step length 1, the probability Pj(c) will be very close to 1 when c≥5 [27]. It implies that 5σz will suppress most of the data in Zjdϵ(x,y). In the present work, we use
(12)Tjdp=51+2ρj,j+12·σjσj+1
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(|d(x,y)|)/0.6745, where d(x,y)∈HHH1 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 d-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 {Y1,Y2,…,YN} of a sequence Y=(Yk)n∈ℕ of independent d-dimensional real random vectors, a lower bound ρ∈[0,∞) for the minimum signal-to-noise ratio (SNR), and a probability value Q≤1-N/(N-2)2. (2) Compute nmin according to h=1/4N(1-Q) and
(13)nmin=N2-hN,
and set λ=2Γ((d+1)/2)/Γ(d/2). (3) Denote by Y(1),Y(2),…,Y(N) the sequence of observations Y1,Y2,…,YN sorted by increasing norm, and compute ξ(ρ)=cosh-1(exp(ρ2/2))=(1/2)ρ+(1/ρ)log(1+1-exp(-ρ2)) if d=1 and ξ(ρ)=(1/ρ)I0-1(exp(ρ2/2)) if d=2, where I0 is the zero-th order modified Bessel function of the first kind. (4) Search a smallest integer n in {nmin,…,N} such that
(14)∥Y(n)∥≤(M{∥Y1∥,∥Y2∥,…,∥YN∥}*(n)λ)ξ(ρ)<∥Y(n+1)∥,
where
(15)M{∥Y1∥,∥Y2∥,…,∥YN∥}*(n)={1n∑k=1n‍∥Y(k)∥ifn≠00ifn=0,
if there exists a such integer n; then set n*=n. Otherwise, set n*=nmin. (5) Compute the estimate σ{Y1,Y2,…,YN}* of the noise standard deviation by using
(16)σ{Y1,Y2,…,YN}*=M{∥Y1∥,∥Y2∥,…,∥YN∥}*(n*)λ.


## 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]
(5)|W2jf(x)|≤K(2j)α,
where j is the decomposition scale, K 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
(6)Pjf(x)=∏i=k1k2‍Wj+if(x),
where k1 and k2 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 k1=0 and let k2=1; then, the DWT multiscale products are
(7)Pjf(x)=Wjf(x)·Wj+1f(x).Obviously, the number of the multiscale products varies with the type of transform. For UDT, the multiscale products have2n1+n2-1 components
(8)Pjdf(x,y)=Wjdf(x,y)·Wj+1df(x,y),
where d 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(Mi,2,1), for i=1,2,…,N, to the logarithmically transformed noisy SAR image f(x,y) up to J scales.(2)
For each coset,
(i)
multiply the directionlet coefficients at adjacent two scales to obtain the multiscale productsPjdf(x,y)=Wjdf(x,y)·Wj+1df(x,y), j=1,2,…,J, d=LHL,LHH,HHL,LLH,HLL,HLH,HHH (L stands for lowpass filtering and H stands for highpass filtering);(ii)
compute the thresholdTjdp and apply it to Pjdf(x,y) to identify W^jdf(x,y) by
(9)W^jdf(x,y)={Wjdf(x,y)-TjdpWjdf(x,y),Pjdf(x,y)≥Tjdp0,Pjdf(x,y)<Tjdp.(3)
Do the inverse S-AWT usingW^jdf(x,y) to reconstruct the f^i(x,y); then perform an exponential transform and obtain the denoised images I^i(x,y).(4)
Compute the final despeckled imageI^(x,y) using I^(x,y)=1/N∑i=1N‍I^i(x,y).It can be seen that the performance of despeckling algorithm depends crucially on the threshold valueTjdp. Thus, the main task of the rest is how to estimate Tjdp 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. LetWjdϵ(x,y) and Wj+1dϵ(x,y) be zero-mean jointly Gaussian distribution with covariance matrix
(10)C=[σj2ρj,j+1σjσj+1ρj,j+1σjσj+1σj+12],
where ρj,j+1 is the correlation coefficient between Wjdϵ(x,y) and Wj+1dϵ(x,y). For the coefficients of directionlets, ρj,j+1 is calculated by ρj,j+1=∫‍∫‍ψj(x,y)·ψj+1(x,y)dxdy/∬(ψj(x,y))2dxdy·∬(ψj+1(x,y))2dxdy, where ψj(x,y) and ψj+1(x,y) are the directionlets used at j and j+1 levels, respectively. The product Zjdϵ(x,y)=Wjdϵ(x,y)·Wj+1dϵ(x,y) has the following PDF [29]:
(11)p(z)=1πσjσj+1exp(ρj,j+1xσjσj+1(1-ρj,j+12))×K0(|z|σjσj+1(1-ρj,j+12)),
where K0 is the zero-th order modified Bessel function of the second kind.The standard deviation ofZjdϵ(x,y) is σz=1+2ρj,j+12·σjσj+1. By computing the values of probability Pj(c)=P{z≤c·σz}, where constant c varies from 1 to 5 by step length 1, the probability Pj(c) will be very close to 1 when c≥5 [27]. It implies that 5σz will suppress most of the data in Zjdϵ(x,y). In the present work, we use
(12)Tjdp=51+2ρj,j+12·σjσj+1
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(|d(x,y)|)/0.6745, where d(x,y)∈HHH1 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 d-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 {Y1,Y2,…,YN} of a sequence Y=(Yk)n∈ℕ of independent d-dimensional real random vectors, a lower bound ρ∈[0,∞) for the minimum signal-to-noise ratio (SNR), and a probability value Q≤1-N/(N-2)2. (2) Compute nmin according to h=1/4N(1-Q) and
(13)nmin=N2-hN,
and set λ=2Γ((d+1)/2)/Γ(d/2). (3) Denote by Y(1),Y(2),…,Y(N) the sequence of observations Y1,Y2,…,YN sorted by increasing norm, and compute ξ(ρ)=cosh-1(exp(ρ2/2))=(1/2)ρ+(1/ρ)log(1+1-exp(-ρ2)) if d=1 and ξ(ρ)=(1/ρ)I0-1(exp(ρ2/2)) if d=2, where I0 is the zero-th order modified Bessel function of the first kind. (4) Search a smallest integer n in {nmin,…,N} such that
(14)∥Y(n)∥≤(M{∥Y1∥,∥Y2∥,…,∥YN∥}*(n)λ)ξ(ρ)<∥Y(n+1)∥,
where
(15)M{∥Y1∥,∥Y2∥,…,∥YN∥}*(n)={1n∑k=1n‍∥Y(k)∥ifn≠00ifn=0,
if there exists a such integer n; then set n*=n. Otherwise, set n*=nmin. (5) Compute the estimate σ{Y1,Y2,…,YN}* of the noise standard deviation by using
(16)σ{Y1,Y2,…,YN}*=M{∥Y1∥,∥Y2∥,…,∥YN∥}*(n*)λ.


## 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 L=3,5, and 9. To simulate the speckle, we adopt the following two-step approach: (1) choose M=1 and generate a standard normal distribution; then determine a one-to-one correspondences between m and L by combining (3) and L=(M/σ)2, where M and σ are the mean and the standard deviation of the log-normal distribution, respectively; (2) then, for a specific L, we use the calculated m 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’ Symmlet4 was employed and the image was decomposed into four resolution levels. Considering the computational complexity, a set of four directionlet transforms S-AWT(Mi,2,1)  (i=1,2) has been used in the proposed method. The corresponding generator matrices are M1=[1001], M2=[11-11], M3=[0110], and M4=[-1111]. 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 Q 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 (S/MSE) ratio to measure the quality of noise suppression. The S/MSE is defined as [10]
(17)SMSEdB=10log10(∑i=1NSi2∑i=1N(S^i-Si)2),
where N is the image size and S and S^ 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)β=Γ(ΔS-ΔS¯,ΔS^-ΔS^¯)Γ(ΔS-ΔS¯,ΔS-ΔS¯)·Γ(ΔS^-ΔS^¯,ΔS^-ΔS^¯),
where ΔS and ΔS^ are the highpass-filtered versions of the original image S and the denoised image S^, respectively, obtained with a window size of 3×3 standard approximation of the Laplacian operator. The overline operator represents the mean value, and Γ(S1,S2)=∑i=1N‍S1i·S2i. The correlation measure, β, should be close to unity for an optimal effect of edge preservation.Figure2(a) illustrates the original test image Syn-SAR. The values of S/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 S/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 (L=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 S/MSE and β measures.Table 1
S/MSE and β values obtained by all despeckling methods applied on the Syn-SAR image.
Method
ENL = 3
ENL = 5
ENL = 9
S/MSE (dB)
β
S/MSE (dB)
β
S/MSE (dB)
β
Speckled
9.6235
0.1140
11.9799
0.1748
16.8858
0.2046
WHT
14.9322
0.1562
15.7528
0.1819
16.9868
0.2283
WMP
15.9533
0.2310
16.9628
0.2808
18.5386
0.3967
Cauchy
15.4109
0.2282
16.9491
0.2963
19.4615
0.4750
Proposed
16.5515
0.2422
17.7322
0.3101
19.7430
0.4861Results obtained by despeckling theSyn-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.
(a)
(b)
(c)
(d)
(e)
(f)


### 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 Q 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 S/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 (I/I^) [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) NoisySuburban 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)(a) NoisyMountain 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)To assess the quality of despeckled images, we show the resultant images in Figures3(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.
Method
WHT
WMP
Cauchy
Proposed
Suburban
ENL
99.2117
103.0979
101.8755
106.5053
Mean
1.0385
1.0364
1.0381
1.0261
Mountain
ENL
45.3968
58.2193
60.2685
60.7680
Mean
1.0323
1.0285
1.0225
1.0204


## 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 L=3,5, and 9. To simulate the speckle, we adopt the following two-step approach: (1) choose M=1 and generate a standard normal distribution; then determine a one-to-one correspondences between m and L by combining (3) and L=(M/σ)2, where M and σ are the mean and the standard deviation of the log-normal distribution, respectively; (2) then, for a specific L, we use the calculated m 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’ Symmlet4 was employed and the image was decomposed into four resolution levels. Considering the computational complexity, a set of four directionlet transforms S-AWT(Mi,2,1)  (i=1,2) has been used in the proposed method. The corresponding generator matrices are M1=[1001], M2=[11-11], M3=[0110], and M4=[-1111]. 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 Q 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 (S/MSE) ratio to measure the quality of noise suppression. The S/MSE is defined as [10]
(17)SMSEdB=10log10(∑i=1NSi2∑i=1N(S^i-Si)2),
where N is the image size and S and S^ 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)β=Γ(ΔS-ΔS¯,ΔS^-ΔS^¯)Γ(ΔS-ΔS¯,ΔS-ΔS¯)·Γ(ΔS^-ΔS^¯,ΔS^-ΔS^¯),
where ΔS and ΔS^ are the highpass-filtered versions of the original image S and the denoised image S^, respectively, obtained with a window size of 3×3 standard approximation of the Laplacian operator. The overline operator represents the mean value, and Γ(S1,S2)=∑i=1N‍S1i·S2i. The correlation measure, β, should be close to unity for an optimal effect of edge preservation.Figure2(a) illustrates the original test image Syn-SAR. The values of S/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 S/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 (L=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 S/MSE and β measures.Table 1
S/MSE and β values obtained by all despeckling methods applied on the Syn-SAR image.
Method
ENL = 3
ENL = 5
ENL = 9
S/MSE (dB)
β
S/MSE (dB)
β
S/MSE (dB)
β
Speckled
9.6235
0.1140
11.9799
0.1748
16.8858
0.2046
WHT
14.9322
0.1562
15.7528
0.1819
16.9868
0.2283
WMP
15.9533
0.2310
16.9628
0.2808
18.5386
0.3967
Cauchy
15.4109
0.2282
16.9491
0.2963
19.4615
0.4750
Proposed
16.5515
0.2422
17.7322
0.3101
19.7430
0.4861Results obtained by despeckling theSyn-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.
(a)
(b)
(c)
(d)
(e)
(f)


## 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 Q 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 S/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 (I/I^) [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) NoisySuburban 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)(a) NoisyMountain 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)To assess the quality of despeckled images, we show the resultant images in Figures3(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.
Method
WHT
WMP
Cauchy
Proposed
Suburban
ENL
99.2117
103.0979
101.8755
106.5053
Mean
1.0385
1.0364
1.0381
1.0261
Mountain
ENL
45.3968
58.2193
60.2685
60.7680
Mean
1.0323
1.0285
1.0225
1.0204


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



---

*Source: 197159-2013-12-02.xml*