16 Information Extraction and Despeckling of SAR Images with Second Generation of Wavelet Transform

Synthetic Aperture Radar (SAR) technology is mainly used to obtain high-resolution images of ground areas in resolutions even less than meter. SAR is even capable of imaging a wide area of terrain and from two and more images it is possible to reconstruct a 3D digital elevation model of ground terrain. Good thing about SAR is an all whether operation and possibility to capture images under various inclination angles. Because digital images are usually corrupted by noise that arises from an imaging device, there is always a need for a good filtering algorithm to remove all disturbances, thus enabling more information extraction. The SAR images are corrupted by a noise called speckle, which makes the interpretation of SAR images very difficult. The goal of removing speckles from the SAR image is to represent a noise-free image and preserve all important features of the SAR image, as for example edges, textures, region borders, etc.


Introduction
Synthetic Aperture Radar (SAR) technology is mainly used to obtain high-resolution images of ground areas in resolutions even less than meter.SAR is even capable of imaging a wide area of terrain and from two and more images it is possible to reconstruct a 3D digital elevation model of ground terrain.Good thing about SAR is an all whether operation and possibility to capture images under various inclination angles.Because digital images are usually corrupted by noise that arises from an imaging device, there is always a need for a good filtering algorithm to remove all disturbances, thus enabling more information extraction.The SAR images are corrupted by a noise called speckle, which makes the interpretation of SAR images very difficult.The goal of removing speckles from the SAR image is to represent a noise-free image and preserve all important features of the SAR image, as for example edges, textures, region borders, etc.
Many different techniques for SAR image despeckling have been proposed over the past few years.Speckle is a noise-like characteristic of SAR images and it is a multiplicative nature, if the intensity or amplitude image is observed.The despeckling can be performed in the image or in the frequency domain.The well-known despeckling filters are Lee (Lee, 1980), Kuan (Kuan et al., 1985), and Frost (Frost et al., 1982).Lee and Kuan filters can be considered as an adaptive mean filters, meanwhile the Frost filter can be considered as a mean adaptive weighted filter.The Bayesian filters are based on the Bayesian theorem, which defines a posterior probability by using a prior, likelihood and evidence probability density functions (pdf).The solution for noise-free image is found by a maximum a posteriori (MAP) estimate.The MAP estimate of a noise free image was proposed in (Walessa & Datcu, 2000), where the noise free image was approximated by a Gauss-Markov random field prior and the noise was modeled with Gamma pdf.Model based despeckling and information extraction is one of the promising techniques of SAR image denoising and scene interpretation.The wavelet based despeckling algorithms have been proposed in (Dai et al., 2004), (Argenti & Alparone, 2002), and (Foucher et al., 2001).The second generation wavelets Chirplet (Cui & Wong, 2006), Contourlet (Chuna et al., 2006), Bandelet (Le Pennec & Mallat, Apr 2005) have appeared over the past few years.
First transform we used is so called Bandelet transform (Le Pennec & Mallat, Dec 2005), which further divides wavelet subbands into smaller subbands using a rate distortion optimization that enables removing redundancy in wavelet transformation.Bandelets (Le Pennec & Mallat, Dec 2005) contain anisotropic wavelets which combine redundancy in the geometric flow of an image corresponding to local directions of its grey levels.With this geometric flow wavelet warping represents a vector field with indication of regularity along edges.Bandelet decomposition is constructed in much the same way as wavelet with use of dyadic squares containing information about bandelet coefficients (parameterized geometric flow) and segmentation (Le Pennec & Mallat, Apr 2005).These squares summarize geometry by local clustering of similar directional vectors.A Bandelet transform can be viewed as an adaptive wavelet basis transform, which is warped according to local direction.
Bandelet transform is therefore capable to separate two different surface areas with different curvatures, which are then decomposed into optimal estimations of regularity direction (Le Pennec & Mallat, Apr 2005).The geometry itself is obtained with regularity flow estimation.Fig. 1 shows an example of directions acquired with bandelet transform.The computational complexity of this transform is much higher as in the case of the classical dyadic decomposition.

Fig. 1. Directions obtained by bandelet transform
The contourlet transform (Do & Vetterli, 2005) is organized a little bit different, because this transform is directly constructed in a discrete space.Thus, contourlet does not need to be transformed from continuous time-space domain.In order to capture as much as possible directional information a 2D directional filter bank is used in contourlet transform (Do & Vetterli, 2005).Directional filter bank is represented with k-binary tree which decomposes original image into 2k bands.These directional filter banks have a flaw mainly because they are designed to capture a direction, which is mainly done in high frequency spectrum of the input image, therefore low frequencies are obstructed.Low frequencies can easily penetrate into several different directional subbands, thus corrupting the transformation subbands.To solve this problem a multiscale decomposition is created with directional decomposition with the help of Laplacian pyramid as a low frequency filter.Laplacian pyramid throughput is a band-pass image which is then led to directional filter bank where a directionality of an image is captured.This scheme can be further applied on a coarser image and thus an iterative scheme can be achieved.It can be concluded that applying iterative contourlet scheme derives to directional subbands in a presence of multiple different scales.
For the despeckling of TerraSAR-X (Wikipedia, 2011) images we used a model based approach, which is supported by first order Bayesian inference.After applying transforms to images a general Gaussian distribution appears in wavelet domain.In this wavelet domain we get subbands different in scales and frequency.The subbands in the wavelet domain have Gaussian distribution and therefore the general Gaussian model is used for a prior density function (pdf).The likelihood pdf is modeled using Gaussian pdf in both, bandelet and contourlet transforms.The despeckling using contourlet (Li et al., 2006) and bandelet (Sveinsson et al., 2008) transforms showed superior despeckling results for SAR images comparing with the wavelet based methods.The model based despeckling mainly depends on chosen models.The image and noise models in the wavelet domain are well defined with presented results in (Argenti & Alparone, 2002), (Gleich & Datcu, 2006) and usually noise-free image is computed using maximum a posteriori estimate.
The despeckling methods were tested using synthetic and real TerraSAR-X data, which were captured in the high resolution spotlight mode.The experimental results showed that the best despeckling method for synthetic images is bandelet transform, because contourlet transform produces artifacts in the homogenous areas.The ratio images between original and despeckled images were examined in order to show estimation of speckle noise, edge and texture preservation using bandelet and contourlet transform.The contourlet transform produces artifacts in form of lines in both homogenous areas and edges.

Second generation wavelets
In this section a comparison between bandelets and contourlets is presented.Bandelets and contourlets are presented in great detail, including subbands creation and filter decomposition.These two denoising schemes are a foundation of later proposed model, which builds a denoising scheme on top of these two schemes yielding better denoising results.

Bandelets
Bandelets (Le Pennec & Mallat, Apr 2005), (Le Pennec & Mallat, Dec 2005) belong to a second generation of wavelet transforms and are composed of anisotropic wavelets, which are in fact a combination of geometric flow of an image corresponding to local directions of its gray levels.This geometric flow represents a regularity of a vector field along edges contained in the image.Typical example of this geometric flow can be seen on Fig. 1, where it can be observed that all directions are aligned to object's edges at the boundary of two different areas.
Edges inside an image are often hard to determine.First generation of bandelet transform uses the vector field (Le Pennec & Mallat, 2001), which determines image regularities and irregularities.Therefore bandelet coefficients represent geometric flow defined by polynomial function.This geometric flow consists of directions of variations in image grey levels, where linear geometric flow is preferred.Bandelet transform image is divided into regions with corresponding vector fields, which describes directions of regularity inside a predefined neighborhood.
If the image intensity is uniformly regular in the neighborhood of a point, then this direction is not uniquely defined, and some form of global regularity is therefore imposed on the flow to specify it uniquely.
In literature it has been proven that the first generation of bandelets has minimum distortion for images whose edges correspond to geometric regularity.However, the first generation of bandelets is composed in continuous space, thus not being able to represent a multiresolution of the geometric regularity.Thus, the second generation of bandelets (Le Pennec & Mallat, Apr 2005) was introduced, which is an orthogonal multiscale transform constructed directly in discrete domain.The bandelet transform first creates a composition of smaller images representing subbands, and then uses fast subband-filtering algorithms.For applications including speckle-noise removal, the geometric flow is optimized in a way that bandelet transform produces minimum distortion in reconstructed images.The decomposition on a bandelet basis is computed using a wavelet filter bank followed by adaptive geometric orthogonal filters, which require O((log 2 N) 3 ) operations.
The key parameters in bandelet transform are: the estimation of basis shapes, the partition of images, and the optimization of geometric flows (Yang et al., 2007).To represent image with as little as possible information, the complex edges must be divided into simpler smaller shapes so that linear geometric flows can represent them sufficiently.The image is commonly divided into smaller square regions that are being divided until there is only one contour inside a square region.It must be noted that the geometric regularity should be discrete, so dyadic decomposition by successive subdivisions of square regions into four smaller sub-squares of twice smaller width can be made.There is a defined maximum and minimum block size (Le Pennec & Mallat, Apr 2005), where the first division produces blocks of maximum size, while later iterations divide those blocks up until minimum size is reached.This partition result can be viewed as a quadtree, where each block is represented by its corresponding leaf in a tree.At each scale the resulting geometry is multiscale and calculated by a procedure that minimizes the Lagrangian cost function.

Implementation
The bandelet transform first computes the 2D wavelet transform of the input original image (Peyré & Mallat, Apr 2005).This transform is based on orthogonal or biorthogonal filter banks and results in four smaller images (children) containing low-and high-frequency components.By selecting a dyadic square and recursively splitting input wavelet image four new sub-squares are created.Further on geometric flow parameterization is performed in each of these sub-squares in every possible direction.Let us assume that each of these squares has a width of k pixels then the number of potential directions d is a little less than 2k 2 .The sampling location is then projected along potential direction and afterwards sorting the resulting 1D points is performed from left to right direction.These points define 1D discrete signal f d (Le Pennec & Mallat, Apr 2005) which is later on transformed with 1D discrete wavelet transform.For a given user defined threshold T, the bandelet transform has to find the best available direction, which in fact produces the less approximation error.Best geometry is obtained by choosing best direction d that minimizes the Lagrangian where f dR is the recovered signal from quantized coefficients acquired by inverse 1D wavelet transform, R G is the number of bits needed to code geometric parameter d, R B is the number of bits needed to code the quantized coefficients and = 3/28 (Le Pennec & Mallat, Apr 2005).
When there are gathered all approximations over each individual dyadic square, the quadtree can be constructed.The algorithm starts with the smallest squares that represent a leaf in quadtree and initialize the cumulative Lagrangian of the sub-tree.

Contourlets
Contourlet transform (Do & Vetterli, 2005) is also classified as a second generation wavelet transform for which a Fourier transform is not needed anymore.Main advantage of second generation wavelet transform over the first generation is a true discrete 2D transformation, which is able to capture geometry of an image, but the first generation wavelet transform does not perform very well on edge regions.This transformation therefore results in adaptive multi-resolution and directional image expansion using contour segments.
Best performance of wavelets is achieved in 1D case which is for example only one row of a 2D picture.Because pictures are not simply stacks of rows, discontinuities evolve along smooth regions.2D wavelet transform thus captures edge points, but on the other hand the throughput on smooth regions is not quite as good anymore.Moreover wavelet transform can only capture a fraction of image directionality, which is clearly seen in Fig. 3 where wavelet transform needs a lot more subdivisions and information then a contourlet transform.In order to capture as much as possible directional information a new type of filter bank has to be constructed.Thus a 2D directional filter bank (Bamberger & Smith, 1992) is used in Contourlet transform.Directional filter bank is represented by k-binary tree, which decomposes original image into 2 k subbands as represented in Fig. 4. The algorithm based on contourlet transform uses a simpler version of directional filter bank, where the first part is constructed from two-channel quincunx filter bank (Vetterli, 1984), while the second part is sampling and reordering operator.With this composition a frequency partition is achieved and also a perfect reconstruction is obtained.As shown in Fig. 6 one can obtain different 2D spectrum decompositions with appropriate combinations of aforementioned building blocks.Thus a k-level binary tree directional filter bank can be viewed as 2 k parallel channel filter bank with equivalent filters and its sampling matrices as shown in Fig. 5  These directional filter banks have a flaw mainly because they are designed to capture directions, which is mainly done in high frequency spectrum of the input image and thus low frequencies are obstructed.As Fig. 5 shows frequency partition a low frequencies can easily penetrate into several different directional subbands and therefore corrupt the transformation subbands.It is therefore wise to combine multiscale decomposition with directional decomposition with the help of Laplacian pyramid as a low frequency filter.Laplacian pyramid throughput is a bandpass image which is then led to directional filter bank where a directionality of the image is captured.This scheme can be further applied on a coarser image and thus an iterative scheme can be achieved.It can be concluded that applying iterative contourlet scheme derives to directional subbands in presence of multiple different scales, which is depicted in Fig. 7.

Bayesian inference incorporated into second generation wavelet transform
The first level of Bayesian inference is given by pyx px pxy py where y represents a noisy image, x represents a noise-reduced image, the 's are the model parameters, p(y|x, ) denotes the conditional pdf called likelihood, p(x| ) denotes prior pdf, and p(y| ) represents evidence pdf.In Eq. ( 2), the evidence pdf does not play a role in the maximization over x, and therefore, the MAP estimator can be written by where the likelihood and prior pdfs should be defined.The MAP estimator is an optimal estimator and minimizes the given cost function.The speckle noise in SAR images is modeled as multiplicative noise, i.e. y = x • z, where z represents pure speckle noise.A multiplicative speckle noise can also be modeled using an additive signal-dependent model, as proposed in (Argenti & Alparone, 2002) where n is a nonstationary signal-dependent additive noise equal to x(z -1).
Models describing texture parameters are widely used in SAR image despeckling (Walessa & Datcu, 2000).Let us model the image as generalized Gauss-Markov random fields (GGMRF) given by

. First order cliques
Let σ x and s define the GGMRF with a neighborhood set s .The MRF model characterizes the spatial statistical dependence of 2D data by a symmetric set called neighboring set.The expression Σ r (x s+r + x s-r ) in Eq. ( 4) represents the sum of all the distinct cliques of neighboring pixel at a specific subband level.A clique is defined as a subset of sites neighboring the observed pixel, where every pair of sites is neighbors of each other.In this double site, cliques are used.A sum is performed over horizontal and vertical neighboring www.intechopen.compixels, for the first model order of the MRF.The neighbor set for a first model order is defined as ζ = {(0, 1), (0, −1), (1, 0), (−1, 0)}, and can be seen in Fig. 8.Moreover, the neighbor set for a second model order is defined as ζ = {(0, 1), (0, −1), (1, 0), (−1, 0), (1, 1), (−1, −1), (1, −1), (−1, 1)}.The MRF model is defined for the symmetric neighbor set; therefore, if r s , then -r s , and ζ is defined as ζ = (r: r s ) ∪ (−r: r s ).The parameter in (4) represents the shape parameter of the GGMRF, Γ(•) represents the Gamma function, and is given by A likelihood pdf is given by a Gaussian distribution where σ n 2 is a noise variance.
The noise variance σ n 2 can be estimated using the results presented in (Argenti & Alparone, 2002), and is given by where x = E[x], and E[x] denotes a mathematical expectation.C x 2 is given by The normalized standard deviation of the noisy wavelet coefficient is given by C Wy = σ Wy / y , where σ Wy is a standard deviation calculated within the wavelet domain, and y is the mean value calculated in the spatial domain.C F denotes the normalized standard deviation of the speckle noise.The parameter ψ l is defined as a product of the coefficients from high-pass (g k ) and low-pass (h k ) filter used at the l-th level of the wavelet decomposition.If the wavelet coefficients of a diagonal detail are of interest, then the parameter ψ l is given by Moreover, if the wavelet coefficients in the horizontal and vertical details are of interest, then the parameter ψ l is given by    The parameter L represents the number of looks of the original SAR image.However, its value is unknown, thus an approximation has to be done, which is L = 2 /σ 2 .The noise variance is then given by where y = E[y].Noise-reduced variance can be computed using the results presented in the paper (Argenti & Alparone, 2002).Thus, noise-reduced the variance is given by Where x 2 is the mean value calculated within the wavelet domain over a predefined window size.
The MAP estimate using the GGMRF primarily defined in ( 4) and the likelihood defined in ( 6) is given by The evidence maximization algorithm is used in order to find the best model's parameters ( , ).The analytical solution for the integral over the posterior p(y|x)p(x| ) does not exists; therefore, the evidence is approximated.The multidimensional pdf is approximated by the multivariate Gaussian distribution with Hessian matrix H centered on the maximum of the a posteriori distribution (Walessa & Datcu, 2000), (Sivia, 1996).The integral over a posterior pdf consists of mutually independent random variables; therefore, a conditional pdf can be rewritten as a product of their components where ∆x = x − x and Hessian matrix H is a square matrix of the second-order partial derivatives of a univariate function The MAP estimate is computed using a numerical method.The texture parameters of the GGMRF model are estimated using the Minimum Mean Square Error (MMSE) estimation technique, and therefore given by a linear model as A logarithmic form can be introduced to simplify Eq. ( 16) as where h ii are the diagonal elements of the Hessian matrix H, which has dimensions of N × N, and N represents the dimension of moving window.Another approximation is then made This approximation is possible, because all off-main diagonal elements represent covariances, and these are sparsely set matrixes that are close to zero, therefore those elements can be neglected.This assumption is in accordance with the statistical independence in Eq. ( 16).Only main-diagonal elements are needed for the Hessian matrix H, which are defined by

Outline of the proposed algorithm
1.The proposed despeckling algorithm transforms SAR images using bandelet (Le Pennec & Mallat, Apr 2005) or contourlet (da Cunha et al., 2006) transform.The number of decompositions depends on the size of the image.The number of levels l is chosen in such a way that the size of the approximation subband is larger than or equal to 64 × 64 pixels (minimum size).2. The model parameters and are estimated inside a window with a size of N × N pixels.In all experiments, a window with 7 × 7 pixels was used.3. The noise and signal variances are estimated using ( 13) and ( 14), respectively.4. The parameter is estimated using the MMSE defined in (18). 5.The shape parameter is changed within the interval [0.5 … 2.5] with a step size of 0.1.6.The noise-reduced coefficients are estimated using the MAP estimate (15) for each value .Each time, the texture parameters are estimated using the MAP estimate obtained in the previous step.
7. The MAP estimate is used for the evidence estimation (21).8.The best MAP estimate x is accepted where the evidence has maximum value.9.The algorithm proceeds to the next pixel.

Synthetic SAR images
The synthetic SAR image, shown in Fig. 11, is composed of four different areas and with added four-look multiplicative speckle noise.The SAR image size shown in Fig. 11  The contourlet transform is constructed using eight directions at the first level of decomposition.The last two levels are chosen to be dyadic, but this is not a requirement.The despeckled images obtained using the bandelet and contourlet transform are shown in Fig. 11 b) and c), respectively.Moreover, the despeckled image obtained using the MBD (Walessa & Datcu, 2000) is shown in Fig. 11 d).
And now with the MAP incorporated into the bandelet and contourlet transform.1, the objective measurements are presented for the denoising of image shown in Fig. 12. Objective measurements include the mean-square error (MSE), the equivalent number of looks (ENL), the mean value of the despeckled image, the ENL of speckle noise (ENL(y/ x )), and the mean value of speckle noise y/ x .The ENL of the image is given by 2 /σ 2 .The best MSE results are from bandelet transform in combination with Bayesian inference, thus having better results than those obtained from the contourlet transforms.A drawback of the contourlet transform is that it produces contours in the reconstructed image, which affects a MSE value.All wavelet-based methods well preserve the mean of the despeckled images.On the other hand, the MBD method well estimates the speckle noise, but it overblurs the reconstructed image, yielding a worse MSE value.
Figs. 13 a)-c) show the ratio images between the original and the reconstructed mosaic images obtained with bandelet, contourlet, and MDB.From these ratio images we can conclude that edges are well preserved, and that the speckle noise in the homogeneous areas is well estimated (i.e.removed) using the MBD method and second generation wavelets, as reported in Table 1.The despeckling within the bandelet and dyadic wavelet domain are able to remove speckles around the strong scatterers, while the contourlet transform produces artifacts in this configuration.Higher image values are difficult to despeckle, because of the nature of the contourlet transform.Therefore, the noise is still present in those areas of the reconstructed image.However, the bandelet transform is overall computationally more demanding than contourlet transform (around 5.6 times), yet the despeckling of each contourlet subband takes about 4.5 times longer than with bandelet transform.Therefore, these methods are also computationally comparable.
To extract texture information from the denoised TerraSAR-X images we have used General Gauss-Markov Random Fields (GGMRF) as a prior pdf (Gleich & Datcu, 2007).As a prior pdf a first order model was used with cliques defined as Gauss-Markov Random Fields and shown in Figs. 9 and 10.Cliques were used to estimate central pixels for both transforms created in a 7 × 7 window which is moving throughout the whole picture.This was applied on transform's first approximation and its corresponding subbands.The texture parameters are iteratively estimated until second order Bayesian inference is increasing, which is used for finding the best model (Gleich & Datcu, 2007).The results of this method can be seen in Fig. 18, where the classification parameters for K-means algorithm were 5 classes and 7 iterations.

Conclusion
This book chapter has presented the proposed methods for despeckling a synthetic and real SAR images using second-generation wavelets.The Bayesian approach in incorporated into second generation wavelets using the wavelet domain.The prior and likelihood pdfs are modeled using GGMRF and Gaussian distribution.The second order Bayesian inference is used to better estimate model parameters and to find the best values possible.The evidence has been simplified and approximated using the Hessian approach.The experimental results have shown that the despeckling of real SAR images using second-generation wavelets is comparable with the dyadic wavelet-based despeckling algorithm (Gleich & Datcu, 2006).Moreover, information extracted using the contourlet domain gives good results using synthetic as well as real SAR data.Unfortunately, the contourlet-based despeckling introduces lines, which are consequences of cutting low-frequency components in the subband decomposition, which can be corrected by introducing a new filter or by post-processing step.

Fig. 2 .
Fig. 2.An example of image denoising using a bandelet transforms.a) Original image, b) Denoised image using a bandelet transforms, and c) Dyadic squares tree Fig. 3. a) Subdivision comparison between wavelet, and b) Contourlet transform

Fig. 7 .
Fig. 7. Construction of contourlet filter bank variable z of the speckle noise is normalized (i.e.E[z] = 1), the parameter C F for intensity images is given by 1

Fig. 10 .
Fig. 10.Neighborhood cliques' organization for contourlet transforms Fig. 11.a) Original speckled image, b) Despeckled image obtained with the original bandelet denoising scheme, c) Despeckled image obtained with the original contourlet denoising scheme, and d) Despeckled image obtained with the MBD denoising technique Fig. 12. a) Original speckled image, b) Despeckled image obtained with the bandeleted denoising scheme, c) Despeckled image obtained with the contourlet denoising scheme, and d) Despeckled image obtained with the MBD denoising technique

Fig. 13 .Fig. 16 .Fig. 17 .
Fig. 13.Ratio images y/ x .a) Ratio image obtained with the bandelet-based despeckling, b) Ratio image obtained with the countourlet-based despeckling, and c) Ratio image obtained with the MBD method The efficiency of the texture separation regarding the proposed method is demonstrated on four Brodatz textures, which are presented in a single mosaic composition and shown in Fig. 14.The textures are corrupted with a four-look speckle noise.The estimated parameter 2 obtained from bandelet and contourlet transforms is shown in Fig. 14 b) and

Fig. 19 .
Fig. 18.Comparison on information extraction.A) Original TerraSAR-X© image, b) Classified image on bandelet transform subbands, and c) Classified image on contourlet transform subbands Texture parameters obtained during the despeckling procedure of the SAR image shown in Fig. 15 with bandelet, contourlet, and MBD method are shown in Fig. 19 b)-d).The algorithm used for classification into four different classes is the K-means algorithm.Fig.19 a) is an indication of K-means algorithm applied to original image scene, where no textures can be identified as no processing was applied.The texture parameters obtained with both proposed algorithms clearly separate between homogeneous and heterogeneous areas.The contourlet transform compared to bandelet transform better separates the homogeneous and heterogeneous areas.From images it can be concluded, that contourlet transform is able to separate more heterogeneous areas from homogeneous ones.As a comparison, the MBD