Undecimated wavelet based speckle reduction for SAR images

Pattern Recognition Letters 26 (2005) 793–800 www.elsevier.com/locate/patrec

Undecimated wavelet based speckle reduction for SAR images D. Gnanadurai a

a,*

, V. Sadasivam

b

Electronics and Communication Engineering Department, Mepco Schlenk Engineering College, Sivakasi, Tamil Nadu 626005, India b Department of Computer Science and Engineering, Manonmanium Sundaranar University, Tirunelveli, Tamil Nadu, India Received 12 October 2003; received in revised form 9 June 2004

Abstract Synthetic Aperture Radar (SAR) images are inherently affected by multiplicative speckle noise, which is due to the coherent nature of scattering phenomena. This paper proposes a novel filtering method for removing such speckle noise from Synthetic Aperture Radar image, that combines the Stationary Wavelet Transform (SWT) with a mean based smoothing operation, which adapts to variation in both signal and the noise. Experimental results on several test images by using the proposed method show that, the proposed method yields significantly superior image quality and better Peak Signal to Noise Ratio (PSNR). In order to illustrate the effectiveness of the proposed algorithm this method is compared with other wavelet based thresholding techniques and conventional statistical filters.
2004 Elsevier B.V. All rights reserved. Keywords: Stationary Wavelet Transform; Filter banks; Bi-orthogonal wavelet; Synthetic Aperture Radar; Speckle noise

1. Introduction Synthetic Aperture Radar (SAR) images are becoming more widely used in many applications such as high-resolution remote sensing for mapping, surface surveillance, search-and-rescue, mine detection, and Automatic Target Detection (ATR). SAR uses microwave radiation to illuminate the earthÕs surface and therefore overcomes the problem associated with conventional visual

*

Corresponding author. E-mail address: [email protected] (D. Gnanadurai).

remote sensing imagery. For example, SAR is not affected by cloud cover and variations in solar illumination. The coherent microwave illumination, however, generates a multiplicative speckle noise that corrupts the images. The major issue on SAR imagery is that basic texture are generally affected by multiplicative speckle noise (e.g., Goodman, 1976). Moreover the presence of speckle damages the radiometric resolution and it affects the task of human interpretation and scene analysis. Thus, it appears sensible to reduce speckle in SAR images, provided that the edges and textural information are not lost. Methods used previously to reduce speckle noise in SAR

0167-8655/$ - see front matter
2004 Elsevier B.V. All rights reserved. doi:10.1016/j.patrec.2004.09.034

794

D. Gnanadurai, V. Sadasivam / Pattern Recognition Letters 26 (2005) 793–800

images are the local statistical methods (e.g., Lee, 1981; Durand and Gimonet, 1987), sigma filter (e.g., Lee, 1983, 1986; Lee and Jurkevich, 1990), median filter (e.g., Lim and Jae, 1990) and homomorphic filtering (e.g., Franceschetti and Pascazio, 1995; Arsenault and Levesque, 1984). In imaging system removal of noise without blurring the edges of image is very critical. Typically, noise is characterized by high spatial frequencies in an image, Fourier based methods usually try to suppress high frequency components, which also tends to affect sharpness of edges. Since the wavelet transform provides good localization in both spatial and spectral domains, low pass filtering is inherent to this transform. Many wavelet based thresholding techniques like hard and soft thresholding were developed by Donoho (1995a,b) and Donoho and Johnstone (1992, 1994, 1995) and all these works were applied for removing White Gaussian Noise only. We describe here an efficient thresholding technique for denoising the speckle, by decomposing the corrupted images into many channels by using Stationary Wavelet Transform (e.g., Cohen et al., 1992). This paper is organized as follows. Section 2 provides a brief introduction to the SAR image and speckle noise. In Section 3, the theory of Stationary Wavelet Transform (undecimated wavelet transform) and its filter bank structure are briefly reviewed. Next in Section 4 the proposed thresholding technique and its shrinkage function are explained. Section 5 presents the performance analysis of the proposed method in denoising the speckle noise for standard image as well as SAR images and it is compared with the standard statistical filters and other wavelet based thresholding techniques. Finally concluding remarks are given in Section 6.

2. SAR image and speckle noise Active radar sensing is often a prime source of inventory information about remote and cloudcovered areas of the world. Due to its high penetration power, Synthetic Aperture Radar acquires high-resolution images in almost all atmospheric

conditions. However, the automatic interpretation of SAR images is often extremely difficult due to speckle noise. Appearing as a random granular pattern, speckle seriously degrades the image quality and hampers the interpretation of image content. Speckle noise in SAR images arises as a consequence of the coherent illumination used by radar. Within each ground resolution cell a large number of elementary scatterers return the radar wave towards the sensor. The backscattered waves being coherent and having different phases undergo a constructive or a destructive interference in a random manner. For a surface that is rough on the scale of the radar wavelength, the number of elementary scatterers within a resolution cell is large enough to ensure the statistical independence in phase and amplitude of the elementary backscattered waves. Theoretically the speckle noise can be modeled by the relation J ¼I þnI

ð1Þ

where n is uniformly distributed random noise with mean 0 and variance v, I is the original Image Matrix and J is the Speckle Noise image.

3. Stationary Wavelet Transform Wavelets are functions generated from one single function W by dilations and translations. The basic idea of the wavelet transform is to represent any arbitrary function as a superposition of wavelets. Any such superposition decomposes the given function into different scale levels, where each level is further decomposed with a resolution adapted to that level. Generally there are two types of wavelet decomposition namely Discrete Wavelet Transform (DWT) (e.g., Daubechies, 1992; Antonini et al., 1992) and Stationary Wavelet Transform (SWT) (e.g., Nanson and Silverman, 1995; Guo, 1995). The DWT is shift variant because, the transform coefficients behave unpredictably under shifts of input signal, a problem that has been treated by introducing large amount of redundancy into the transform to make it shift invariant. In SWT a transform-domain redundancy of 4.0 has been introduced to make it shift invariant. Here

D. Gnanadurai, V. Sadasivam / Pattern Recognition Letters 26 (2005) 793–800

in this work SWT is used, because of its following merits over DWT: (i) Even though the SWT requires more calculation and calls for bigger memory, it enables better denoising quality and better edge detecting capacity. (ii) Since the number of coefficients is halved in subsequent resolution levels of DWT, due to down sampling (decimation) process, it is usually suitable for implementing for discrete signals or images whose size is power of 2. But SWT can be implemented for any arbitrary size of images, since the down sampling process is not applied so as to keep the number of coefficients is same in all the resolution levels. For this reason, the SWT can also be called as undecimated wavelet transform. This SWT is normally finds its application in image texture analysis and enhancement techniques. Due to the Decomposition of an image using the Stationary Wavelet Transform the original image is transformed into four pieces which can be labeled as LL1, LH1, HL1, HH1 as in the schematic depicted in Fig. 1. The LL piece comes from low pass filtering in both horizontal and vertical directions and it is most like the original picture and so is called the course or approximated sub-band. The remaining pieces LH1, HL1 and HH1 are called as detailed components. The HL1 comes from

Columns LD

LL1

HD

LH1

LD

HL1

HD

HH1

Row s LD

X

795

Columns LL1

LR

Rows LR

LH1

HR

HL1

LR

X

HR HH1

HR

Fig. 2. SWT reconstruction filter bank structure.

low pass filtering in the vertical direction and high pass filtering in the horizontal direction and so has the label HL1. The visible detail in this sub-image HL1, such as edges, have an overall vertical orientation since their alignment is perpendicular to the direction of the high pass filtering and they are called vertical details. Similarly, the LH1 is constructed from low pass filtering in the horizontal direction and high pass filtering in the vertical direction and contains horizontal details. Since the HH1 component is the result of high pass filtering in both vertical and horizontal directions, it contains the diagonal details only. The filters LD and HD shown in Fig. 1 are onedimensional Low Pass Filter (LPF) and High Pass Filter (HPF) respectively for image decomposition. To obtain the next level of decomposition, sub-band LL1 alone is further decomposed. This process continues until some final scale is reached. In this method if an image of size 256 · 256 is decomposed up to 3 levels, 10 sub-images of the same 256 · 256 size are obtained. The decomposed image can be reconstructed using a reconstruction filter as shown in Fig. 2. Here, the filters LR and HR represents low pass and high pass reconstruction filters respectively.

HD

4. Proposed method

Fig. 1. SWT decomposition filter bank structure.

The block diagram shown in Fig. 3 depicts the various steps involved in the proposed denoising

796

D. Gnanadurai, V. Sadasivam / Pattern Recognition Letters 26 (2005) 793–800

Finest Bands

Noisy Image

Stationary Wavelet Decomposition

Mean Based Smoothing Operation

Stationary Wavelet Reconstruction

De-noised Image

Average Filtering Course Band Fig. 3. Block diagram of proposed work.

method. This filtering method for removing speckle noise combines the Stationary Wavelet Transform with a mean-based smoothing operation. In this work the noisy image is a first decomposed up to five levels using bi-orthogonal wavelet (bior 2.2) to obtain 16 sub-bands. Then in each subband, except the lowest band (coarse coefficients), mean-based smoothing operation is applied by using proposed shrinkage function, which is defined mathematically for an image I(x, y) as
Iðx; yÞ ¼

Iðx; yÞ M

if jIðx; yÞj < M if jIðx; yÞj > M

ð2Þ

where M is the mean value of row of the pixel I(x, y) and x, y—pixel coordinates. The above representation can be explained as follows. First mean value (M) of each row is found out. Then all the wavelet coefficients in that row are compared with this mean value, and if the absolute value of a coefficient in the row exceeds the mean value then it is replaced by the mean value itself and if the value of wavelet coefficient is less than the mean value then the corresponding value is left unchanged. This operation is smoothening the sudden changes of values of coefficients, which may be due to random noise. The same operation is applied to all rows of the sub-bands. The graphical representation of the shrinkage function given in Eq. (2) is depicted in Fig. 4. Since the lowest band (LL band) consists of courser pixel value of image, average filter is applied to this band to remove noise. After this smoothing operation all the sub-bands are combined by inverse Stationary Wavelet Transform to reconstruct the denoised image.

Fig. 4. The shrinkage function of proposed thresholding technique.

5. Experimental results In this section, the simulation results, obtained by processing several standard test images using our proposed method are presented and the results of the proposed approach are compared with other current state-of-the-art speckle filtering methods. This method is also compared with DWT based soft and hard thresholding methods (DonohoÕs universal thresholding technique). In order to quantify the achieved performance improvement, two different measures are computed based on original and denoised image. For quantitative evaluation, an extensively used measure is Peak Signal to Noise Ratio (PSNR) (e.g., Amir Said, 1994) which is denoted as   2552 PSNR ¼ 10log10 dB ð3Þ MSE

Table 1 PSNR values of different denoised test images Noise variance

Noise image

Median filter

Wiener filter

Hard threshold (DWT)

Soft threshold (DWT)

Proposed SWT method

Lena Barbara Seed Pepper

0.04

19.80 21.60 19.77 19.61

24.81 26.12 25.12 24.20

25.86 26.84 26.04 24.99

26.23 26.53 27.28 24.95

27.61 27.10 29.24 26.15

28.41 26.76 29.30 25.81

Lena Barbara Seed Pepper

0.08

16.93 18.62 16.91 16.97

22.13 23.63 22.23 21.70

23.15 23.95 23.24 22.25

24.30 24.16 24.90 22.90

26.01 25.66 27.03 24.60

27.04 26.16 27.58 24.85

Lena Barbara Seed Pepper

0.12

15.30 16.95 15.26 15.51

20.43 22.10 20.58 20.13

21.54 22.38 21.68 20.65

23.07 22.84 23.71 21.71

24.84 24.80 25.48 23.50

25.98 25.62 26.21 24.01

Lena Barbara Seed Pepper

0.16

14.15 15.76 14.11 14.46

19.26 20.93 19.39 19.14

20.42 21.23 20.63 19.57

22.14 21.96 22.79 20.73

23.85 24.15 24.50 22.64

25.10 25.16 25.12 23.22

Table 2 S/M values of different denoised test images Noise variance

Noise image

Median filter

Wiener filter

Hard threshold (DWT)

Soft threshold (DWT)

Proposed SWT method

Lena Barbara Seed Pepper

0.04

0.4388 0.5338 0.3718 0.6866

0.4005 0.4898 0.3298 0.6623

0.3923 0.4869 0.3207 0.6513

0.3896 0.4916 0.3167 0.6509

0.3818 0.4742 0.3099 0.6410

0.3772 0.4656 0.3049 0.6338

Lena Barbara Seed Pepper

0.08

0.4809 0.5779 0.4170 0.7145

0.4164 0.5049 0.3460 0.6728

0.3998 0.4996 0.3275 0.6548

0.3892 0.4960 0.3180 0.6510

0.3808 0.4749 0.3067 0.6370

0.3756 0.4683 0.3015 0.6302

Lena Barbara Seed Pepper

0.12

0.5176 0.6145 0.4557 0.7386

0.4285 0.5178 0.3605 0.6836

0.4047 0.5093 0.3328 0.6586

0.3908 0.5007 0.3159 0.6490

0.3789 0.4755 0.3041 0.6340

0.3737 0.4679 0.2972 0.6264

Lena Barbara Seed Pepper

0.16

0.5510 0.6493 0.4901 0.7618

0.4423 0.5325 0.3733 0.6929

0.4112 0.5203 0.3364 0.6621

0.3893 0.5042 0.3158 0.6485

0.3770 0.4744 0.3003 0.6309

0.3705 0.4677 0.2916 0.6232

797

Image (512 · 512)

D. Gnanadurai, V. Sadasivam / Pattern Recognition Letters 26 (2005) 793–800

Image (512 · 512)

798

D. Gnanadurai, V. Sadasivam / Pattern Recognition Letters 26 (2005) 793–800

Fig. 5. De-noising of Lena image corrupted by Speckle noise of variance 0.08: (a) noise image (16.93 dB); (b) proposed method (27.04 dB).

where the Mean Square Error (MSE) is given by MSE ¼

K 1 X 2 ðX  P Þ K i¼1

ð4Þ

where X is the original image, P is the denoised image and K is the image size. The performance of the proposed method is checked with many standard images like Lena, Barbara, Pepper and Seed. Speckle noise with different variance is applied to these images artificially and the corrupted image is denoised by the proposed method as well as the statistical filters like median and wiener filters and wavelet based thresholding techniques like soft thresholding and

Table 3 S/M values of different denoised SAR images Denoising methods

SAR 1

SAR 2

Noisy image Median filter Wiener filter Hard thresholding Soft thresholding Proposed method

0.5619 0.4687 0.4735 0.5131 0.4924 0.4121

0.6201 0.5068 0.4860 0.5051 0.4710 0.4358

SAR 1—ÔHorse TrackÕ image (Fig. 6). SAR 2—ÔMaricopa agricultural centerÕ image (Fig. 7).

hard thresholding. The obtained results (PSNR in dB) are tabulated in Table 1. From Table 1, it is

Fig. 6. Processing of SAR image of ÔHorse TrackÕ (near Albuquerque, NM): (a) original image (S/M = 0.5619); (b) denoised image (S/M = 0.4121).

D. Gnanadurai, V. Sadasivam / Pattern Recognition Letters 26 (2005) 793–800

799

Fig. 7. Processing of SAR image of the view of Maricopa agricultural center (near Phonix, AZ): (a) original image (S/M = 0.6201); (b) denoised image (S/M = 0.4358).

obvious the proposed SWT method yields higher PSNR than other DWT based thresholding methods and statistical filters. Also, in order to quantify the speckle reduction performance the standard deviation (S) to mean (M) ratio (S/M) (e.g., Alin et al., 2003) is computed for denoised images. This quantity is a measure of image speckle in homogeneous regions. Standard deviation to mean ratio (S/M) is also calculated for standard images, which measures the speckle reduction performance. These values are tabulated in Table 2 and from this table it is obvious the proposed method yields minimum value, which indicates high speckle reduction. Fig. 5(b) shows the denoised image of noisy image of variance 0.08 which PSNR is 16.93 dB and for this case the proposed method yields 27.04 dB. It is clear, from this enhanced Lena image (Fig. 5(b)), the noise reduction capability by using the proposed SWT method is remarkable. In order to prove the merit of proposed method in practical applications, noisy SAR images are chosen which is affected by speckle noise during acquisition and the proposed algorithm is applied without adding artificial noise. The S/M values of denoised SAR images are tabulated in Table 3. The S/M values of the denoised images by proposed method are minimum, which means high speckle suppression. The first SAR image shown in Fig. 6(a) depicts a ÔHorse trackÕ image (near Albuquerque, New Mexico) affected by speckle phenomena which S/M is equal to 0.5619. The denoised image by the proposed method is depicted in Fig. 6(b) which S/M value is reduced to 0.4121.

The one more SAR image shown in Fig. 7(a) depicts the view of Maricopa agricultural center (near Phonix, AZ), affected by speckle phenomena which S/M is 0.6201. The denoised image by the proposed method is depicted in Fig. 7(b) which S/M value is reduced to 0.4349.

6. Conclusion This paper presents an effective and low-complexity, Stationary Wavelet Transform based method for speckle reduction. Since the proposed method is based on the analysis of mean value of the wavelet coefficients, it is more adaptive to sub-bands. The experimental results of proposed method are compared with other methods, in order to illustrate the effectiveness of the proposed algorithm. The comparison suggests that the proposed method results are competitive and outperforming with other wavelet based methods and statistical filters. Since this proposed method proves better results in standard images, it can be used for denoising of images obtained from coherent imaging systems as holographic ones, Synthetic Aperture Radar (SAR) and ultrasound imaging systems, which output image is always affected by speckle noise.

References Alin, A., Panagiotius, T., Anastasios, B., 2003. SAR image denoising via Bayesian wavelet shrinkage based on

800

D. Gnanadurai, V. Sadasivam / Pattern Recognition Letters 26 (2005) 793–800

heavy-tailed modeling. IEEE Trans. Geosci. Remote Sensing 41 (8), 1779–1782. Amir Said, 1994. A new and efficient image codec based on set partitioning in hierarchical tress. IEEE Trans. Circuit System Video Technol. 6, 48. Antonini, M., Barlaud, M., Mathieu, P., Daubechies, I., 1992. Image coding using wavelet transform. IEEE Trans. Image Process. 1 (2), 205–220. Arsenault, H.H., Levesque, M., 1984. Combined homomorphic and local-statistics processing for restoration of images degraded by signal-dependent noise. Appl. Opt. 23 (6), 845– 850. Cohen, A., Daubechies, I., Feauveau, J., 1992. Bi-orthogonal bases of compactly supported wavelets. Commn. Pure Appl. Math. 45, 485–560. Daubechies, I., 1992. Ten Lectures on Wavelets. SIAM, Philadelphia, PA. Donoho D.L., Johnstone I.M., 1994. Ideal spatial adaptation by wavelet shrinkage. Biometrika, pp. 425–455. Donoho, D.L., Johnstone, I.M., 1992. Minimax estimation via wavelet shrinkage. Technical Report, Statistics, Stanford. Donoho, D.L., 1995a. De-noising by soft-thresholding. IEEE Trans. Inform. Theory 41, 613–627. Donoho, D.L., 1995b. Nonlinear solution of linear inverse problems by Wavelet Decomposition. Appl. Comput. Harm. Anal. 2, 101–126. Donoho, D.L., Johnstone, I.M., 1995. Adapting to un-known smoothness via wavelet shrinkage. J. Amer. Stat. Assoc. 90 (432), 1200–1224.

Durand, J.M., Gimonet, B.J., 1987. SAR data filtering for classification. IEEE Trans. Geosci. Remote Sensing 25 (5), 629–637. Franceschetti, G., Pascazio, V., 1995. Iterative homomorphic technique for speckle reduction in synthetic-aperture radar imaging. J. Opt. Soc. Amer. 12 (4), 686–694. Goodman, J.W., 1976. Some fundamental properties of speckle. J. Amer. 66, 1145–1150. Guo, H., 1995. Theory and applications of shift invariant, timevarying and undecimated wavelet transform. Master Thesis. Rice University, Houston, TX. Lee, J.S., 1981. Speckle analysis and smoothing of synthetic aperture radar images. Comput. Graphics Image Graphics 17, 24–32. Lee, J.S., 1983. A simple speckle smoothing algorithm for synthetic aperture radar images. IEEE Trans. Systems Man Cybernet. 13 (1), 85–89. Lee, J.S., 1986. Speckle suppression and analysis for synthetic aperture radar images. Opt. Engrg. 25 (5), 636–643. Lee, J.S., Jurkevich, I., 1990. Coastline detection and tracing in SAR images. IEEE Trans. Geosci. Remote Sensing 28 (4), 662–668. Lim, Jae, S., 1990. Two-Dimensional Signal and Image Processing. Prentice-Hall, Englewood Cliffs, NJ, pp. 469– 476. Nanson, G.P., Silverman, B.W., 1995. Stationary wavelet transform and some statistical applications. Technical Report, Department of Mathematics, University of Bristol, Bristol, UK.