**2.2 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.

Information Extraction and Despeckling of

2 1 *E <sup>l</sup>* 

different scales, which is depicted in Fig. 7.

Fig. 7. Construction of contourlet filter bank

sampling matrix

*E*0

*E*1

SAR Images with Second Generation of Wavelet Transform 379

0 *y*

<sup>0</sup> *S*

<sup>1</sup> *S*

2 1 *S <sup>l</sup>* 

*x x*ˆ

1 *y*

2 1 *y <sup>l</sup>* 

*x x*ˆ

1 *y*

Fig. 6. 2D partition of spectrum using quincunx filter banks with fan filters. Darker shades represent the ideal frequency supports of each filter. Denotion Q represents a quincunx

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

0 *y*

<sup>0</sup> *S*

<sup>1</sup> *S*

2 1 *S <sup>l</sup>* 

Fig. 5. The multichannel view of a *k*-level tree-structured directional filter bank

*D*0

*D*1

2 1 *D <sup>l</sup>* 

Fig. 3. a) Subdivision comparison between wavelet, and b) 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 (Do & Vetterli, 2005). In Fig. 5 *D* denotes an equivalent directional filter.

Fig. 4. Frequency partitioning where *k* = 3 and there are 23 = 8 real wedge-shaped frequency bands. Subbands 0-3 correspond to the mostly horizontal directions, while subbands 4-7 correspond to mostly vertical directions

a) b)

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 (Do

> 2

> > 2

1

2

5 6

Fig. 4. Frequency partitioning where *k* = 3 and there are 23 = 8 real wedge-shaped frequency bands. Subbands 0-3 correspond to the mostly horizontal directions, while subbands 4-7

1

 ,

4

5

7

0

3

1

Fig. 3. a) Subdivision comparison between wavelet, and b) Contourlet transform

& Vetterli, 2005). In Fig. 5 *D* denotes an equivalent directional filter.

0

3

4

6

7

 ,

correspond to mostly vertical directions

Fig. 5. The multichannel view of a *k*-level tree-structured directional filter bank

Fig. 6. 2D partition of spectrum using quincunx filter banks with fan filters. Darker shades represent the ideal frequency supports of each filter. Denotion Q represents a quincunx sampling matrix

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.

Fig. 7. Construction of contourlet filter bank

#### **3. Bayesian inference incorporated into second generation wavelet transform**

The first level of Bayesian inference is given by

$$p\left(\mathbf{x}|y,\theta\right) = \frac{p\left(y|\mathbf{x},\theta\right)p\left(\mathbf{x}|\theta\right)}{p\left(y|\theta\right)}\tag{2}$$

Information Extraction and Despeckling of

SAR Images with Second Generation of Wavelet Transform 381

pixels, 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

<sup>1</sup> <sup>3</sup> , <sup>1</sup> *x x*

<sup>2</sup>

The noise variance σn2 can be estimated using the results presented in (Argenti & Alparone,

2 22 2 1

*C* 

*l F*

 2 21 2 2 *<sup>l</sup> lk k k k*

*h g*

Moreover, if the wavelet coefficients in the horizontal and vertical details are of interest,

Since the random variable *z* of the speckle noise is normalized (i.e. E[z] = 1), the parameter

*l kk k k*

*h g*

2 1 2 2 *<sup>l</sup>*

<sup>2</sup> 1 *Wy l F*

2 2

*C C*

exp 2 2

1

*y x py x* 

 

*s s*

 *n lx F x C C*

where *μx* = E[x], and E[x] denotes a mathematical expectation. Cx2 is given by

2

*x*

coefficients of a diagonal detail are of interest, then the parameter *ψl* is given by

The normalized standard deviation of the noisy wavelet coefficient is given by *CWy* = *σWy*/*μy*, where *σWy* is a standard deviation calculated within the wavelet domain, and *μy* is the mean value calculated in the spatial domain. *CF* denotes the normalized standard deviation of the speckle noise. The parameter *ψl* is defined as a product of the coefficients from high-pass (*gk*) and low-pass (*hk*) filter used at the *l*-th level of the wavelet decomposition. If the wavelet

*C*

A likelihood pdf is given by a Gaussian distribution

where σn2 is a noise variance.

then the parameter *ψl* is given by

*CF* for intensity images is given by

2002), and is given by

2 2

*n n*

(5)

(6)

(7)

(8)

(9)

(10)

1 *C L <sup>F</sup>* (11)

*s s*

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

$$
\hat{\mathfrak{x}}(y) = \arg\max\_{\mathfrak{x}} p(y|\mathbf{x}, \boldsymbol{\theta}) p(\mathbf{x}|\boldsymbol{\theta}) \tag{3}
$$

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) *y = x · z = x + x(z – 1) = x + n*, 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

$$p\left(\mathbf{x}\_{s}|\boldsymbol{\theta}\right) = \frac{\nu\eta\left(\boldsymbol{\nu}, \sigma\_{\boldsymbol{x}}\right)}{2\Gamma\left(1/\nu\right)} \exp\left[-\left[\eta\left(\boldsymbol{\nu}, \sigma\_{\boldsymbol{x}}\right)\middle|\mathbf{x}\_{s} - \sum\_{r \in \mathcal{L}\_{s}} \theta\_{r}\left(\mathbf{x}\_{s+r} + \mathbf{x}\_{s-r}\right)\right]\right]^{\boldsymbol{\nu}}\tag{4}$$

Fig. 8. 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(xs+r + xs-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 pixels, 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

$$\eta\left(\nu,\sigma\_{\ge}\right) = \sigma\_{\ge}^{-1}\sqrt{\frac{\Gamma\left(3/\nu\right)}{\Gamma\left(1/\nu\right)}}\tag{5}$$

A likelihood pdf is given by a Gaussian distribution

$$p\left(y\_s | \mathbf{x}\_s\right) = \frac{1}{\sqrt{2\pi\sigma\_n^2}} \exp\left(-\frac{\left(y\_s - \mathbf{x}\_s\right)^2}{2\sigma\_n^2}\right) \tag{6}$$

where σn2 is a noise variance.

380 Advances in Wavelet Theory and Their Applications in Engineering, Physics and Technology

**3. Bayesian inference incorporated into second generation wavelet transform** 

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

> <sup>ˆ</sup> arg max , *<sup>x</sup> x y p yx p x*

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) *y = x · z = x + x(z – 1) = x + n*, where *n* is a non-

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

, exp , 2 1 *<sup>s</sup>*

1

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(xs+r + xs-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

2

*r*

 

*s x s r sr sr*

*p x x xx*

Neighborhood

Cliques

Parameters

*p xy*

the maximization over *x*, and therefore, the MAP estimator can be written by

stationary signal-dependent additive noise equal to *x(z – 1).*

  

*x*

(GGMRF) given by

Fig. 8. First order cliques

, , *p yx p x*

 (2)

(3)

(4)

*p y* 

The first level of Bayesian inference is given by

The noise variance σn2 can be estimated using the results presented in (Argenti & Alparone, 2002), and is given by

$$
\sigma\_n^2 = \varphi\_1 \mu\_x^2 \mathbf{C}\_F^2 \left( \mathbf{1} + \mathbf{C}\_x^2 \right) \tag{7}
$$

where *μx* = E[x], and E[x] denotes a mathematical expectation. Cx2 is given by

$$\mathbf{C}\_x^2 = \frac{\mathbf{C}\_{\text{My}}^2 - \boldsymbol{\nu}\_I \mathbf{C}\_F^2}{\boldsymbol{\nu}\_I \left(\mathbf{1} + \mathbf{C}\_F^2\right)}\tag{8}$$

The normalized standard deviation of the noisy wavelet coefficient is given by *CWy* = *σWy*/*μy*, where *σWy* is a standard deviation calculated within the wavelet domain, and *μy* is the mean value calculated in the spatial domain. *CF* denotes the normalized standard deviation of the speckle noise. The parameter *ψl* is defined as a product of the coefficients from high-pass (*gk*) and low-pass (*hk*) 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

$$\mathcal{W}\_l = \left(\sum\_k h\_k^2\right)^2 \left(\sum\_k g\_k^2\right)^{2(l-1)}\tag{9}$$

Moreover, if the wavelet coefficients in the horizontal and vertical details are of interest, then the parameter *ψl* is given by

$$\boldsymbol{\nu}\_{l} = \left(\sum\_{k} h\_{k}^{2}\right) \left(\sum\_{k} \mathbf{g}\_{k}^{2}\right)^{2l-1} \tag{10}$$

Since the random variable *z* of the speckle noise is normalized (i.e. E[z] = 1), the parameter *CF* for intensity images is given by

$$C\_F = 1/\sqrt{L} \tag{11}$$

while for the amplitude images the parameter *CF* is given by

$$C\_F = \sqrt{(4/\pi - 1)/L} \tag{12}$$

Information Extraction and Despeckling of

weights the clique on different subbands.

θ1

Ap

θ<sup>2</sup> θ<sup>3</sup>

Fig. 9. Neighborhood cliques' organization for bandelet transforms

A logarithmic form can be introduced to simplify Eq. (16) as

1

*i*

*p yx*

*N*

Ve

θ4

technique, and therefore given by a linear model as

SAR Images with Second Generation of Wavelet Transform 383

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

<sup>1</sup> *T T*

(18)

θ6

*GG G X*

where *X* are the observed coefficients inside the window of a size *N* × *N*, and matrix *G* consists of neighboring coefficients attributed around each individual observed coefficient *xs* inside a window of a size *N* × *N*. The organization of a neighborhood for the bandelet and contourlet domain is shown in Fig. 9 and Fig. 10, respectively. Those figures show the parent-child relationships for the bandelet and contourlet transform. Each parameter *θ*

Ho

Di

θ<sup>8</sup> θ<sup>7</sup>

1

*N ii i H h* 

*ii i i i*

(19)

(20)

*h py x px*

<sup>1</sup> log log 2 log log log ˆ ˆ <sup>2</sup>

where *hii* 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

θ5

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

$$
\sigma\_n^2 = \frac{\mathbb{C}\_F^2 \left(\mathbb{Y}\_l \mu\_y^2 + \sigma\_{\mathbb{W}y}^2\right)}{1 + \mathbb{C}\_F^2} \tag{13}
$$

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

$$
\sigma\_{\mathbf{x}}^{2} = \psi\_{1} \mu\_{\mathbf{x}}^{2} \mathbf{C}\_{\mathbf{x}}^{2} \tag{14}
$$

Where *μx2* 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

$$-\nu\eta \left(\nu\_{\prime}\sigma\_{x}\right) \left| \eta \left(\nu\_{\prime}\sigma\_{x}\right) \right| \mathbf{x}\_{s} - \sum\_{r \neq \zeta\_{s}} \theta\_{r} \left(\mathbf{x}\_{s+r} + \mathbf{x}\_{s-r}\right) \right|^{\nu-1} + \frac{y\_{s} - \mathbf{x}\_{s}}{\sigma\_{n}^{2}} = 0 \tag{15}$$

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

$$\begin{aligned} p\left(y|\boldsymbol{\theta}\right) &= \int p\left(y|\mathbf{x}\right) p\left(\mathbf{x}|\boldsymbol{\theta}\right) d\mathbf{x} \\ p\left(y|\boldsymbol{\theta}\right) &\approx \int \prod\_{i=1}^{N} p\left(y\_{i}|\hat{\mathbf{x}}\_{i}\right) p\left(\hat{\mathbf{x}}\_{i}|\boldsymbol{\theta}\right) \exp\left(-\frac{1}{2}\Delta\mathbf{x}^{T}H\Delta\mathbf{x}\right) d\mathbf{x} \\ p\left(y|\boldsymbol{\theta}\right) &\approx \frac{\left(2\pi\right)^{N/2}}{\sqrt{|H|}} \prod\_{i=1}^{N} p\left(y\_{i}|\mathbf{x}\_{i}\right) p\left(\mathbf{x}\_{i}|\boldsymbol{\theta}\right) \end{aligned} \tag{16}$$

where ∆x = x − *x*ˆ and Hessian matrix *H* is a square matrix of the second-order partial derivatives of a univariate function

$$H = -\nabla\nabla\sum\_{i=1}^{N} \log\left(p\left(y\_i|\mathbf{x}\_i\right)p\left(\mathbf{x}\_i|\theta\right)\right) \tag{17}$$

 4 1 *C L <sup>F</sup>* 

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

where *μy* = E[y]. Noise-reduced variance can be computed using the results presented in the

2 22

Where *μx2* is the mean value calculated within the wavelet domain over a predefined

The MAP estimate using the GGMRF primarily defined in (4) and the likelihood defined in

<sup>2</sup> , , 0

<sup>1</sup> ˆ ˆ exp <sup>2</sup>

*p y py x px x H x dx*

where ∆x = x − *x*ˆ and Hessian matrix *H* is a square matrix of the second-order partial

log *<sup>N</sup> ii i <sup>i</sup> H py x px*

 

*ii i*

 

<sup>1</sup>

*r n y x x xx*

*s s x x s r sr sr*

2

*C*

*n*

paper (Argenti & Alparone, 2002). Thus, noise-reduced the variance is given by

 

*p y p y x p x dx*

*N*

1 2

2

*i*

*i*

*N N*

*p y py x px <sup>H</sup>*

1

*ii i*

 

*s*

 

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

 

2 22

*C*

*F*

 

<sup>2</sup> 1 *F l y Wy*

(12)

(13)

*x lx x C* (14)

(15)

1

*T*

(16)

(17)

while for the amplitude images the parameter *CF* is given by

variance is then given by

window size.

(6) is given by

 

rewritten as a product of their components

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

$$\boldsymbol{\theta} = \left(\boldsymbol{G}\boldsymbol{G}^T\right)^{-1} \left(\boldsymbol{G}^T\boldsymbol{X}\right) \tag{18}$$

where *X* are the observed coefficients inside the window of a size *N* × *N*, and matrix *G* consists of neighboring coefficients attributed around each individual observed coefficient *xs* inside a window of a size *N* × *N*. The organization of a neighborhood for the bandelet and contourlet domain is shown in Fig. 9 and Fig. 10, respectively. Those figures show the parent-child relationships for the bandelet and contourlet transform. Each parameter *θ* weights the clique on different subbands.

Fig. 9. Neighborhood cliques' organization for bandelet transforms

A logarithmic form can be introduced to simplify Eq. (16) as

$$\log p\left(y|\mathbf{x}\right) \approx \sum\_{i=1}^{N} \frac{1}{2} \left( \log \left(2\pi\right) - \log \left(h\_{ii}\right) \right) + \log p\left(y\_i|\hat{\mathbf{x}}\_i\right) + \log p\left(\hat{\mathbf{x}}\_i|\boldsymbol{\theta}\right) \tag{19}$$

where *hii* 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

$$\left| H \right| \approx \prod\_{i=1}^{N} h\_{ii} \tag{20}$$

Information Extraction and Despeckling of

**5. Experimental results 5.1 Synthetic SAR images** 

9. The algorithm proceeds to the next pixel.

7. The MAP estimate is used for the evidence estimation (21).

SAR Images with Second Generation of Wavelet Transform 385

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 is 512 × 512 pixels; therefore three levels of decomposition are used for bandelet transform. First let us show the difference between the pure bandelet and contourlet, and the MBD method.

a) b)

c) d) 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

8. The best MAP estimate *x*ˆ is accepted where the evidence has maximum value.

Fig. 10. Neighborhood cliques' organization for contourlet transforms

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

$$h\_{ii} = \sum\_{i=1}^{N} -\upsilon(\upsilon - 1)\eta\left(\upsilon, \sigma\_x\right)^2 \left(\eta\left(\upsilon, \sigma\_x\right) \left|\hat{\mathbf{x}}\_i - \sum\_{i \in \mathcal{L}} \theta\_j\left(\mathbf{x}\_i^{\prime} + \mathbf{x}\_i^{\prime\prime}\right)\right|\right)^{\upsilon - 2} - \sum\_{i=1}^{N} \frac{\hat{\mathbf{x}}\_i}{\sigma\_n^2} \tag{21}$$

#### **4. Outline of the proposed algorithm**

