Stochastic Systems: Modeling, Optimization, and ApplicationsView this Special Issue
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 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.
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”  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 , 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 . 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 . (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 , 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)) . 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 , 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  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
Since the statistical properties of speckle noise have been studied by Goodman , a large number of models were proposed to analyse the SAR images, such as negative exponential distribution, gamma distribution , log-normal distribution , and generalized gamma distribution . In , 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 log-normal random variable  can be generated using where and are the mean and the median values of the distribution, respectively, and is a standard normal random variable. The equivalence between (number of looks) in a speckle image and in (3) has been given in .
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 where , , , and .
According to lattice theory , any cubic lattice can be partitioned into cosets of the lattice , where each coset is determined by the shift vectors , for . 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 . 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 , 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 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 .
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 . The relation between wavelet coefficients and Lipschitz component satisfies  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 and are nonnegative integers.
In , Bao and Zhang have pointed out that it is sufficient to implement the multiplication of two adjacent scales in practice; thus, we let and let ; 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 components 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  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, for , 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 , , (L stands for lowpass filtering and H stands for highpass filtering);(ii)compute the threshold and apply it to to identify by (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 .
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 be zero-mean jointly Gaussian distribution with covariance matrix where is the correlation coefficient between and . For the coefficients of directionlets, is calculated by , where and are the directionlets used at and levels, respectively. The product has the following PDF : where is the zero-th order modified Bessel function of the second kind.
The standard deviation of is . By computing the values of probability , where constant varies from 1 to 5 by step length 1, the probability will be very close to when . It implies that will suppress most of the data in . In the present work, we use 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 , where 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)  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 for the minimum signal-to-noise ratio (SNR), and a probability value . (2) Compute according to and and set . (3) Denote by the sequence of observations sorted by increasing norm, and compute if and if , where is the zero-th order modified Bessel function of the first kind. (4) Search a smallest integer in such that where if there exists a such integer ; then set . Otherwise, set . (5) Compute the estimate of the noise standard deviation by using
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  (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 , and we considered three different levels of simulated speckle noise corresponding to , and . To simulate the speckle, we adopt the following two-step approach: (1) choose 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 to generate the simulated speckle.
In the WHT, WMP, Cauchy, and the proposed algorithm, the orthogonal wavelet of Daubechies’ Symmlet was employed and the image was decomposed into four resolution levels. Considering the computational complexity, a set of four directionlet transforms S-AWT has been used in the proposed method. The corresponding generator matrices are , , , and . Notice that the four transforms are oriented along , , , and , respectively. Moreover, in the process of estimating noise standard deviation, two parameters and involved in DATE are set as and , respectively.
To evaluate the performances of these despeckling approaches, we computed the signal-to-mean squared error () ratio to measure the quality of noise suppression. The is defined as  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  where and are the highpass-filtered versions of the original image and the denoised image , respectively, obtained with a window size of standard approximation of the Laplacian operator. The overline operator represents the mean value, and . 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 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 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 () 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 and measures.
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 are shown in Figures 3(a) and 4(a), respectively. Since the noise-free images are not available, the values of the 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 ()  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.
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.
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.
V. S. Frost, J. A. Stiles, K. S. Shanmugan, and J. C. Holtzman, “A Model for radar images and its application to adaptive digital filtering of multiplicative noise,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 4, no. 2, pp. 157–166, 1982.View at: Publisher Site | Google Scholar
J. S. Lee, “Digital image enhancement and noise filtering by use of local statistics,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 2, no. 2, pp. 165–168, 1980.View at: Google Scholar
L. Gagnon and A. Jouan, “Speckle filtering of SAR images: a comparative study between complex-wavelet-based and standard filters,” in Wavelet Applications in Signal and Image Processing V, vol. 3169 of Proceedings of the SPIE, pp. 80–91, San Diego, Calif, USA, July, 1997.View at: Publisher Site | Google Scholar
F. Argenti, T. Bianchi, A. Lapini, and L. Alparone, “Simplified MAP despeckling based on Laplacian-Gaussian modeling of undecimated wavelet coefficients,” in Proceedings of the 19th European Signal Processing Conference (EUSIPCO '11), pp. 1140–1144, Barcelona, Spain, Aughust, 2011.View at: Publisher Site | Google Scholar
C. Oliver and S. Quegan, Understanding Synthetic Aperture Radar Images, SciTech, Raleigh, NC, USA, 2004.
J. H. Conway and N. J. A. Sloane, Sphere Packings, Lattices and Groups, vol. 290, Springer, New York, NY, USA, 1998.
Y. Meyer, Wavelets and Operators, vol. 37 of Cambridge Studies in Advanced Mathematics, Cambridge University Press, Cambridge, UK, 1993.
H.-Y. Gao, “Wavelet shrinkage denoising using the non-negative garrote,” Journal of Computational and Graphical Statistics, vol. 7, no. 4, pp. 469–488, 1998.View at: Google Scholar
M. K. Simon, Probability Distributions Involving Gaussian Random Variables: A Handbook for Engineers, Scientists and Mathematicians, Springer, New York, NY, USA, 2006.