Multiscale Filtering of SAR Images Using Scale and ... - CiteSeerX
Multiscale Filtering of SAR Images Using Scale and Space Consistency Samuel Foucher Research & Development Dept Computer Research Institute of Montreal Montreal, Canada Abstract—A new approach for speckle reduction in SAR images based on the stationary wavelet transform is proposed. Noisy wavelet coefficients are reduced via a shrinkage function that depends on a statistical modeling of the normalized wavelet coefficient probability density functions. Consistencies between coefficients across scales and space are also reinforced using consistency rules. The approach is particularly robust in cases of correlated speckle noise. Keywords-speckle filtering; SAR; wavelet;
I. INTRODUCTION Most speckle reduction techniques assume the speckle noise as a white random process. Recently, numerous wavelet based methods have been proposed [1][2][3][4][5] where, like most of the mono-scale filters, a white speckle noise is assumed. However, the impact of the speckle second order statistics is expected to be significant as a non uniform speckle power spectral density function will strongly affect the distribution of the signal energy within the various wavelet sub-bands [5]. The proposed method is based on the work of Scharcanski et al. [6] on the filtering of optical images corrupted by Gaussian noise. In this approach, the amplitude of significant wavelet coefficients is modeled probabilistically and a shrinkage function is derived based on the model obtained. Scale consistency is applied and leads to a better preservation of strong coefficients present at successive scales. Identically, Geometric constraint (orientation preservation) aims at strengthening coefficients placed along edges. The present paper is an extension of the authors previous work on speckle filtering [1][5] where a multiscale MAP estimation was proposed (SWT-MMSE). The paper is organized as follows. Section II, considers the impact of speckle correlation on the normalized wavelet coefficient statistics. In Section III, a wavelet coefficient shrinkage function based on a two components mixture model is described. Finally, in section III the proposed approach is compared with other speckle reduction techniques on simulated and real images.
II.
INFLUENCE OF A CORRELATED SPECKLE ON THE WAVELET COEFFICIENTS
A. General Considerations In the following, we assume a stationary wavelet transform (SWT) of a SAR image I composed of 3J high-frequency images
{ A[ j ]I }
j =1,...,J
{Wε[ j ]I }ε=h,v,d j =1,...,J
and
J
low-frequency
images
. An interesting property of the wavelet
[j ] of the transform is that the autocorrelation function (acf) RW
wavelet coefficients W [ j ]I is the convolution of the acf RI of the SAR image I by the acf Rψ[ j ] of the wavelet filter ψ[ j ] : [j] W[j]I = I ∗ ψ[ j ] ⇒ RW = RI ∗ Rψ[ j ]
(1)
Note that the wavelet acf is also a wavelet. Consequently, the speckle noise power for a L looks SAR image within a homogeneous area of mean reflectivity µI which is expressed by the following [j] = varW I
µI2 ( ρ ∗ Rψ[ j ] ) (0, 0) L S
(2)
The linear filtering power gain is directly affected by the interaction of the speckle correlation function ρS and the wavelet acf, hence ∆
S 2,[ jC] = ( ρS ∗ Rψ[ j ] ) (0, 0) =
1 2 γ (ω) Ψ[ j ](ω) d ω (3) 2π ∫ S
where γS is the speckle power spectral density function and Ψ[ j ] is the wavelet Fourier transform. More particularly, for a white speckle ( γS (ω) = 1 ) the above relationship is simply the wavelet power gain [1]: [j ]
[j]
S 2,C = S 2 =
2
∑ ( ψk[ j ] )
(4)
k
For a correlated speckle noise, the impact will be an increased noise level on lower frequency bands (j>1) and a decrease on finest scale (j=1).
This work has been supported in part by the NSERC of Canada (Discovery Grant) and the MDEIE of the “Gouvernement du Québec”.
B. Properties of the normalized wavelet coefficients We propose the following two operators robust to speckle noise: L M I = [j] A I [j ]
∑
S2,[ jε]
ε = h ,v,d
W [ j ]I θ[ j ]I = arctan v [ j ] W I h
(5)
(6)
1, if (x,y) ∈ edge g [ j ](x , y ) = 0, if (x,y) ∈ noise
The amplitude operator M [ j ]I is an average of the three directional wavelet coefficient energies after gain normalization. Additionally, because of the multiplicative nature of the speckle noise, we normalized M [ j ]I using the low pass image AI [ j ] . Consequently, for a white speckle noise, the mean value of M [ j ]I should be approximately 1. The orientation operator θ[ j ]I indicates the energy direction and will be used below in application of geometric constraint. The proposed approach is based on the fact that M [ j ]I behaves differently for noise and edge related wavelet coefficients. [j]
We assume that M I is the root square of a sum of three independent centered random variables which means that the probability density function (pdf) should be close to a scaled χ distribution with 3 degrees of freedom (noted σ[ j ]χ3 ):
[j]
[j ]
p (M I = r | n) =
( ) e 3 2σ Γ ( ) 2 r σ[ j ] [j ]
2
−
( )
1 r 2 σ[ j ]
2
(7)
[j ]
is the scaling factor that should be where σ approximately 1 for a white speckle noise and around S2,[ jC] / S 2[ j ] for a correlated noise. We verify the validity of this model by simulating a 256x256 uniform speckled image (L=1) with and without correlation. The correlated speckle was generated using the Blacknell technique [8].
Figure 1. Observed pdf of M I for an image of white speckle noise (left) and correlated speckle noise (right). The red curve is the proposed model .
(8)
We assume that the observed pdf for a realization r of the random variable M [ j ]I is a mixture of two pdfs: p[ j ](r ) = π[ j ]p[ j ](r | n) + (1 − π[ j ] )p[ j ](r | e)
(9)
Additionally, we assume scaled chi distributions for the two M [ j ]I | n ∼ σn[ j ]χ3 and mixture components: M [ j ]I | e ∼ σe[ j ]χ3 . There are a variety of ways to estimate the
scaling factors σn[ j ] and σe[ j ] . We suggest using either the mean value or the mode position: σˆn[ j ] = arg max { h [ j ](b) } / 2 b
σˆn[ j ]
3 = Γ( ) M [ j ]I /( 2Γ(2)) 2
(10)
where h [ j ] is the histogram of the observed M [ j ]I . Using the mode position is likely to offer a more robust approach for the estimation of σn[ j ] because it is less sensitive to the tail of the histogram. Conversely, the mean based estimator should be more reliable for σe[ j ] . We designed a quick computation scheme based on two assumptions, where: 1) σˆn[ j ] is determined from the mode of the histogram of M [ j ]I ; and 2) σˆe[ j ] should satisfy the following relation derived from (9) 3 Γ( ) M [ j ]I /( 2Γ(2)) = π[ j ]σˆn[ j ] + (1 − π[ j ] )σˆe[ j ] 2 derived from (9). The pseudo-code is provided in Appendix. The shrinkage function is then the posterior probability of having an edge given the observed M [ j ]I value: g [ j ](r ) = p[ j ](e | r ) =
[j]
PROPOSED APPROACH
A. Shrinkage function The ideal shrinkage function retains wavelet coefficients generated by edges or targets and put to zero those generated by noise:
2
(Wε[ j ]I )
III.
(1 − π[ j ] )p[ j ](r | e) (11) π[ j ]p[ j ](r | n) + (1 − π[ j ] )p[ j ](r | e)
This approach is tested on a simulated SAR image shown on Fig. 3. The estimated mixtures and corresponding shrinkage functions are given in Fig. 2. Since the speckle is correlated, the mode position is greater than one and should fluctuate according to the quantity S2,[ jC] / S 2[ j ] . Consequently, this approach is robust to noise correlation because the resulting shrinkage function will adjust accordingly. Point targets will produce very large wavelet coefficients that must be preserved so we impose g [ j ] = 1 if M [ j ]I > 2σn[ j ] L + 2 .
N α[i ]g [ j ],sc (x + i, y ) ∑ i =−N N ∑ α[i ]g [ j ],sc (x , y + i ) i =−N g [ j ],geo (x , y ) = N ∑ α[i ]g [ j ],sc (x + i, y + i) i =−N N ∑ α[i ]g [ j ],sc (x + i, y − i ) i =−N
, if θ [ j ]I = 00 , θ [ j ]I = 900
(13) [j ]
, θ I = 45
0
, θ [ j ]I = 1350
where α[i ] is a weighting coefficient (Gaussian), θ [ j ]I is a Figure 2. Estimated mixture components (left) and corresponding shrinkage function (right) for three wavelet levels. The green thick curves correspond to
p[ j ](r | e)
and
g[ j ]
respectively.
B. Scale and Space Consistency Wavelet coefficients produced by edges and other meaningful structures within the image are persistent across scales, whereas noise related coefficients have their energy spread evenly on all scales. Consequently, meaningful coefficients will be consistent along scales. This observation can be used in order to increase the signal to noise ratio on the noisy fine scales, Scharcanski et al. [6] proposed the following rule g [ j ],sc (x , y ) =
J − j +1 (12) 1 1 1 + … + g [ j ](x , y ) g [ j +1](x , y ) g [J ](x , y )
This combination rule has the following behaviour: •
If one of the shrinkage functions takes a zero value then the resulting combination will be zero. As a result, homogeneous areas are preserved across scales;
•
For edges, values must be consistent across scales in order not to be penalized (i.e. g [ j ],sc = 1 ⇔ g[ j ] = g [ j +1] = ...g [J ] = 1 ). Therefore, the combination rule (12) is a conjunctive fusion rule ( g [ j ],sc ≈ min { g [n ] } ). n≥j
The geometric consistency rule reinforces wavelet coefficients placed along edges and penalizes isolated coefficients. Edges are assumed to be piecewise linear so wavelet coefficients are locally smoothed according to their energy direction (6) as following,
quantized version of θ[ j ]I and 2N+1 is the size of the neighbourhood. Because we expect edge direction to be consistent in two consecutives scales, we apply an additional direction consistency rule [6] only where g [ j ],geo (x , y ) < 1 in order to preserve point targets: g [ j ],dir (x, y ) = g [ j ],geo (x, y ) ⋅ cos(θ[ j ]I (x, y ) − θ[ j +1]I (x, y )) (14)
C. Algorithm Summary The full non parametric filtering approach is composed of the following steps: 1) Compute the stationary wavelet transform in J levels, obtaining 3J high frequency images and 3J low frequency images: {Wε[ j ]I }
j =1,...,J
ε = h ,v,d
, { A[ j ]I }
j =1,...,J
;
2) Apply the following from the coarsest level to the finest level, a.
calculate the wavelet coefficient magnitudes M [ j ]I (5) and orientations θ[ j ]I (6); b. estimate σn[ j ], σe[ j ] and π[ j ] using the procedure given in the appendix; c. calculate the intra-scale shrinkage function g [ j ] (12); d. Apply the scale consistency rule (12), obtaining g [ j ],sc ; e. apply the geometric consistency rule (13), obtaining g [ j ],geo ; f. apply the direction consistency rule (14), obtaining g [ j ],dir ; g. finally, modify the wavelet coefficients by multiplying by g [ j ],dir ; 3) Apply the inverse wavelet transform from the modified wavelet coefficients, obtaining the filtered image.
IV.
TABLE I.
ENL AND EP VALUES FOR EACH FILTERING RESULT. THE LOWER THE EP VALUE, THE BETTER THE EDGE PRESERVATION.
RESULTS
A. Artificial SAR Image The algorithm is tested on an artificial image shown on Figure 3. The correlated speckle is generated using Blacknell’s algorithm [8] with a large correlation coefficient value ( ρ0,1 = ρ1,0 = 0.45 ). After filtering we estimate the Equivalent Number of Looks (ENL) within the image homogeneous areas, the PSNR and an Edge Preservation (EP) indicator [11] (preservation of the original gradient values along edges). Four levels of decomposition are used for the stationary wavelet transform with Biorthogonal-5 filters. For the proposed method, the following parameters are used: ∆α = 5.10−3 , αmin = 0.1 , αmax = 0.95 and N=3.
Filters
ENL
EP
PSNR
g [ j ] (11)
2.61
0.051
32.9
g [ j ],sc (12)
4.81
0.084
35.2
g [ j ],geo (13)
4.75
0.101
35.4
g [ j ],dir (14)
6.21
0.233
34.3
Refined Lee (9x9) [9]
7.38
0.422
29.3
SWT-MMSE
3.36
0.079
31.6
Proposed Method
B. Real SAR Image The filters are applied on a 550 × 400 E-SAR image of Oberpfaffenhofen, Germany (L-band, hh) (see Fig. 5). We estimate the ENL within an homogeneous area (white square). The proposed method using g [ j ],sc gives 17.9 looks, the SWTMMSE achieved 9.3 looks and Refined Lee, 15.9.
Figure 3. Simulated SAR image with correlated speckle noise (left) and ground truth (right).
Figure 4. Top row: propose method with the shrinkage function (left) and
g [ j ],geo
g [ j ],sc
(right). Bottom row: SWT-MMSE (left) and Refined Lee Filter 9x9 (right).
Figure 5. Original Image (Top Left), proposed method with g [ j ],sc (top right), SWT-MMSE (bottom left) and Refined Lee 9x9 (bottom right). Images are displayed in db.
V. CONCLUSION In this paper, we proposed a new approach for the filtering of SAR images based on the stationary wavelet transform. The denoising of wavelet coefficients is achieved by a shrinkage function derived from a two component mixture model of the pdf of the normalized wavelet coefficients. This approach is particularly suitable for the filtering of correlated speckle and does not require any knowledge of the speckle second-order statistics. Different shrinkage functions are also proposed that are reinforcing consistent values and geometry across scales. Filtering results on artificial and real SAR images show a more suitable compromise between noise reduction and detail preservation as compared to standard filters. Noise reduction is stronger to what can be achieved by other wavelet based speckle reduction techniques that are assuming a white speckle noise. Among the four different shrinkage functions proposed, g [ j ],geo produces a speckle reduction similar to the Enhanced Lee filter (9x9) but with superior point target and edge preservation. The additional direction consistency rule (14) is not improving filtering results because point targets are altered. ACKNOWLEDGMENT This work has been supported in part by the NSERC of Canada (Discovery Grant) and the MDEIE of the “Gouvernement du Québec”. The authors wish to thank Elaine Rosenberg and Lisa Hollinger for their linguistic expertise. APPENDIX The following pseudo-code describes the estimation procedure for σn[ j ] and σe[ j ] (see section III.A). Initialization H ← histogram(M [ j ]I ) σˆn[ j ] ← arg max{H }/ 2 3 σ[ j ] ← ( ) M [ j ]I /( 2Γ(2)) 2 m ←1 α ← αmin While α < αmax σe[ j ](m ) ← (σ[ j ] − α(m )σˆn[ j ] )/(1 − α(m )) RMS (m ) ← computeRMS (H , p[ j ](r )) m ← m +1 α(m ) ← α(m − 1) + ∆α End σˆe[ j ] ← σe[ j ](arg min{RMS }) m
[j ]
πˆ
← α(arg min{RMS }) m
REFERENCES [1]
S. Foucher, G. B. Bénié, and J.-M. Boucher, Multiscale MAP Filtering of SAR images, IEEE Trans. on Image Proc., vol. 10, no. 1, pp. 49 –6, Jan. 2001. [2] A. Achim, P. Tsakalides and A. Bezerianos; “SAR Image Denoising via Bayesian Wavelet Shrinkage Based on Heavy-Tailed Modeling”, IEEE Trans. on Geosci. and Remote Sens., vol 41, no 8, pp. 1773 –1784, Aug. 2003. [3] Hua Xie, L. E. Pierce and F. T. Ulaby, “SAR Speckle Reduction using Wavelet Denoising and Markov Random Field Modeling”, IEEE Trans. on Geosci. and Remote Sens., vol. 40, no. 10, pp. 2196 –2212, Oct. 2002. [4] F. Argenti and L. Alparone, “Speckle Removal from SAR Images in the Undecimated Wavelet Domain”, IEEE Trans. on Geosci. and Remote Sens., vol. 40, no. 11, pp. 2363 –2212, Nov. 2002. [5] S. Foucher, J.-M. Boucher and G. B. Bénié. "Wavelet Filtering of Correlated Speckle" In AeroSense 2003 (SPIE #5108). Orlando, 21-25 April 2003. [6] J. Scharcanski; C. R. Jung, and R. T. Clarke, “Adaptive Image Denoising Using Scale and Space Consistency”, IEEE Trans. on Image Proc., vol. 11, no. 9, pp. 1092 –1101, Sept. 2002 [7] A. Papoulis, Probability, “Random Variables, and Stochastic Processes,” McGraw-Hill, 1991. [8] D. Blacknell, “New method for the simulation of correlated k-distributed clutter,” IEE Proceedings, radar, Sonar Navigation 141, p. 45, Feb. 1994. [9] J.-S. Lee, M.R. Grunes, and G. de Grandi, “Polarimetric SAR Speckle Filtering and Its Implication for Classification,” IEEE Trans. on Geosci. and Remote Sens., vol. 37, no. 5, pp. 2363-2373, 1999. [10] S. Foucher, G. Farage, and G. B. Bénié, “SAR Image Filtering Based on the Stationary Contourlet Transform”, IGARSS, Denver, pp. 4021-4024 Jul. 2006. [11] G. Farage, S. Foucher, and G. B. Bénié, “Comparison of PolSAR Speckle Filtering Techniques”, IGARSS, Denver, pp. 1760-1763, Jul. 2006.
I. INTRODUCTION Most speckle reduction techniques assume the speckle noise as a white random process. Recently, numerous wavelet based methods have been proposed [1][2][3][4][5] where, like most of the mono-scale filters, a white speckle noise is assumed. However, the impact of the speckle second order statistics is expected to be significant as a non uniform speckle power spectral density function will strongly affect the distribution of the signal energy within the various wavelet sub-bands [5]. The proposed method is based on the work of Scharcanski et al. [6] on the filtering of optical images corrupted by Gaussian noise. In this approach, the amplitude of significant wavelet coefficients is modeled probabilistically and a shrinkage function is derived based on the model obtained. Scale consistency is applied and leads to a better preservation of strong coefficients present at successive scales. Identically, Geometric constraint (orientation preservation) aims at strengthening coefficients placed along edges. The present paper is an extension of the authors previous work on speckle filtering [1][5] where a multiscale MAP estimation was proposed (SWT-MMSE). The paper is organized as follows. Section II, considers the impact of speckle correlation on the normalized wavelet coefficient statistics. In Section III, a wavelet coefficient shrinkage function based on a two components mixture model is described. Finally, in section III the proposed approach is compared with other speckle reduction techniques on simulated and real images.
II.
INFLUENCE OF A CORRELATED SPECKLE ON THE WAVELET COEFFICIENTS
A. General Considerations In the following, we assume a stationary wavelet transform (SWT) of a SAR image I composed of 3J high-frequency images
{ A[ j ]I }
j =1,...,J
{Wε[ j ]I }ε=h,v,d j =1,...,J
and
J
low-frequency
images
. An interesting property of the wavelet
[j ] of the transform is that the autocorrelation function (acf) RW
wavelet coefficients W [ j ]I is the convolution of the acf RI of the SAR image I by the acf Rψ[ j ] of the wavelet filter ψ[ j ] : [j] W[j]I = I ∗ ψ[ j ] ⇒ RW = RI ∗ Rψ[ j ]
(1)
Note that the wavelet acf is also a wavelet. Consequently, the speckle noise power for a L looks SAR image within a homogeneous area of mean reflectivity µI which is expressed by the following [j] = varW I
µI2 ( ρ ∗ Rψ[ j ] ) (0, 0) L S
(2)
The linear filtering power gain is directly affected by the interaction of the speckle correlation function ρS and the wavelet acf, hence ∆
S 2,[ jC] = ( ρS ∗ Rψ[ j ] ) (0, 0) =
1 2 γ (ω) Ψ[ j ](ω) d ω (3) 2π ∫ S
where γS is the speckle power spectral density function and Ψ[ j ] is the wavelet Fourier transform. More particularly, for a white speckle ( γS (ω) = 1 ) the above relationship is simply the wavelet power gain [1]: [j ]
[j]
S 2,C = S 2 =
2
∑ ( ψk[ j ] )
(4)
k
For a correlated speckle noise, the impact will be an increased noise level on lower frequency bands (j>1) and a decrease on finest scale (j=1).
This work has been supported in part by the NSERC of Canada (Discovery Grant) and the MDEIE of the “Gouvernement du Québec”.
B. Properties of the normalized wavelet coefficients We propose the following two operators robust to speckle noise: L M I = [j] A I [j ]
∑
S2,[ jε]
ε = h ,v,d
W [ j ]I θ[ j ]I = arctan v [ j ] W I h
(5)
(6)
1, if (x,y) ∈ edge g [ j ](x , y ) = 0, if (x,y) ∈ noise
The amplitude operator M [ j ]I is an average of the three directional wavelet coefficient energies after gain normalization. Additionally, because of the multiplicative nature of the speckle noise, we normalized M [ j ]I using the low pass image AI [ j ] . Consequently, for a white speckle noise, the mean value of M [ j ]I should be approximately 1. The orientation operator θ[ j ]I indicates the energy direction and will be used below in application of geometric constraint. The proposed approach is based on the fact that M [ j ]I behaves differently for noise and edge related wavelet coefficients. [j]
We assume that M I is the root square of a sum of three independent centered random variables which means that the probability density function (pdf) should be close to a scaled χ distribution with 3 degrees of freedom (noted σ[ j ]χ3 ):
[j]
[j ]
p (M I = r | n) =
( ) e 3 2σ Γ ( ) 2 r σ[ j ] [j ]
2
−
( )
1 r 2 σ[ j ]
2
(7)
[j ]
is the scaling factor that should be where σ approximately 1 for a white speckle noise and around S2,[ jC] / S 2[ j ] for a correlated noise. We verify the validity of this model by simulating a 256x256 uniform speckled image (L=1) with and without correlation. The correlated speckle was generated using the Blacknell technique [8].
Figure 1. Observed pdf of M I for an image of white speckle noise (left) and correlated speckle noise (right). The red curve is the proposed model .
(8)
We assume that the observed pdf for a realization r of the random variable M [ j ]I is a mixture of two pdfs: p[ j ](r ) = π[ j ]p[ j ](r | n) + (1 − π[ j ] )p[ j ](r | e)
(9)
Additionally, we assume scaled chi distributions for the two M [ j ]I | n ∼ σn[ j ]χ3 and mixture components: M [ j ]I | e ∼ σe[ j ]χ3 . There are a variety of ways to estimate the
scaling factors σn[ j ] and σe[ j ] . We suggest using either the mean value or the mode position: σˆn[ j ] = arg max { h [ j ](b) } / 2 b
σˆn[ j ]
3 = Γ( ) M [ j ]I /( 2Γ(2)) 2
(10)
where h [ j ] is the histogram of the observed M [ j ]I . Using the mode position is likely to offer a more robust approach for the estimation of σn[ j ] because it is less sensitive to the tail of the histogram. Conversely, the mean based estimator should be more reliable for σe[ j ] . We designed a quick computation scheme based on two assumptions, where: 1) σˆn[ j ] is determined from the mode of the histogram of M [ j ]I ; and 2) σˆe[ j ] should satisfy the following relation derived from (9) 3 Γ( ) M [ j ]I /( 2Γ(2)) = π[ j ]σˆn[ j ] + (1 − π[ j ] )σˆe[ j ] 2 derived from (9). The pseudo-code is provided in Appendix. The shrinkage function is then the posterior probability of having an edge given the observed M [ j ]I value: g [ j ](r ) = p[ j ](e | r ) =
[j]
PROPOSED APPROACH
A. Shrinkage function The ideal shrinkage function retains wavelet coefficients generated by edges or targets and put to zero those generated by noise:
2
(Wε[ j ]I )
III.
(1 − π[ j ] )p[ j ](r | e) (11) π[ j ]p[ j ](r | n) + (1 − π[ j ] )p[ j ](r | e)
This approach is tested on a simulated SAR image shown on Fig. 3. The estimated mixtures and corresponding shrinkage functions are given in Fig. 2. Since the speckle is correlated, the mode position is greater than one and should fluctuate according to the quantity S2,[ jC] / S 2[ j ] . Consequently, this approach is robust to noise correlation because the resulting shrinkage function will adjust accordingly. Point targets will produce very large wavelet coefficients that must be preserved so we impose g [ j ] = 1 if M [ j ]I > 2σn[ j ] L + 2 .
N α[i ]g [ j ],sc (x + i, y ) ∑ i =−N N ∑ α[i ]g [ j ],sc (x , y + i ) i =−N g [ j ],geo (x , y ) = N ∑ α[i ]g [ j ],sc (x + i, y + i) i =−N N ∑ α[i ]g [ j ],sc (x + i, y − i ) i =−N
, if θ [ j ]I = 00 , θ [ j ]I = 900
(13) [j ]
, θ I = 45
0
, θ [ j ]I = 1350
where α[i ] is a weighting coefficient (Gaussian), θ [ j ]I is a Figure 2. Estimated mixture components (left) and corresponding shrinkage function (right) for three wavelet levels. The green thick curves correspond to
p[ j ](r | e)
and
g[ j ]
respectively.
B. Scale and Space Consistency Wavelet coefficients produced by edges and other meaningful structures within the image are persistent across scales, whereas noise related coefficients have their energy spread evenly on all scales. Consequently, meaningful coefficients will be consistent along scales. This observation can be used in order to increase the signal to noise ratio on the noisy fine scales, Scharcanski et al. [6] proposed the following rule g [ j ],sc (x , y ) =
J − j +1 (12) 1 1 1 + … + g [ j ](x , y ) g [ j +1](x , y ) g [J ](x , y )
This combination rule has the following behaviour: •
If one of the shrinkage functions takes a zero value then the resulting combination will be zero. As a result, homogeneous areas are preserved across scales;
•
For edges, values must be consistent across scales in order not to be penalized (i.e. g [ j ],sc = 1 ⇔ g[ j ] = g [ j +1] = ...g [J ] = 1 ). Therefore, the combination rule (12) is a conjunctive fusion rule ( g [ j ],sc ≈ min { g [n ] } ). n≥j
The geometric consistency rule reinforces wavelet coefficients placed along edges and penalizes isolated coefficients. Edges are assumed to be piecewise linear so wavelet coefficients are locally smoothed according to their energy direction (6) as following,
quantized version of θ[ j ]I and 2N+1 is the size of the neighbourhood. Because we expect edge direction to be consistent in two consecutives scales, we apply an additional direction consistency rule [6] only where g [ j ],geo (x , y ) < 1 in order to preserve point targets: g [ j ],dir (x, y ) = g [ j ],geo (x, y ) ⋅ cos(θ[ j ]I (x, y ) − θ[ j +1]I (x, y )) (14)
C. Algorithm Summary The full non parametric filtering approach is composed of the following steps: 1) Compute the stationary wavelet transform in J levels, obtaining 3J high frequency images and 3J low frequency images: {Wε[ j ]I }
j =1,...,J
ε = h ,v,d
, { A[ j ]I }
j =1,...,J
;
2) Apply the following from the coarsest level to the finest level, a.
calculate the wavelet coefficient magnitudes M [ j ]I (5) and orientations θ[ j ]I (6); b. estimate σn[ j ], σe[ j ] and π[ j ] using the procedure given in the appendix; c. calculate the intra-scale shrinkage function g [ j ] (12); d. Apply the scale consistency rule (12), obtaining g [ j ],sc ; e. apply the geometric consistency rule (13), obtaining g [ j ],geo ; f. apply the direction consistency rule (14), obtaining g [ j ],dir ; g. finally, modify the wavelet coefficients by multiplying by g [ j ],dir ; 3) Apply the inverse wavelet transform from the modified wavelet coefficients, obtaining the filtered image.
IV.
TABLE I.
ENL AND EP VALUES FOR EACH FILTERING RESULT. THE LOWER THE EP VALUE, THE BETTER THE EDGE PRESERVATION.
RESULTS
A. Artificial SAR Image The algorithm is tested on an artificial image shown on Figure 3. The correlated speckle is generated using Blacknell’s algorithm [8] with a large correlation coefficient value ( ρ0,1 = ρ1,0 = 0.45 ). After filtering we estimate the Equivalent Number of Looks (ENL) within the image homogeneous areas, the PSNR and an Edge Preservation (EP) indicator [11] (preservation of the original gradient values along edges). Four levels of decomposition are used for the stationary wavelet transform with Biorthogonal-5 filters. For the proposed method, the following parameters are used: ∆α = 5.10−3 , αmin = 0.1 , αmax = 0.95 and N=3.
Filters
ENL
EP
PSNR
g [ j ] (11)
2.61
0.051
32.9
g [ j ],sc (12)
4.81
0.084
35.2
g [ j ],geo (13)
4.75
0.101
35.4
g [ j ],dir (14)
6.21
0.233
34.3
Refined Lee (9x9) [9]
7.38
0.422
29.3
SWT-MMSE
3.36
0.079
31.6
Proposed Method
B. Real SAR Image The filters are applied on a 550 × 400 E-SAR image of Oberpfaffenhofen, Germany (L-band, hh) (see Fig. 5). We estimate the ENL within an homogeneous area (white square). The proposed method using g [ j ],sc gives 17.9 looks, the SWTMMSE achieved 9.3 looks and Refined Lee, 15.9.
Figure 3. Simulated SAR image with correlated speckle noise (left) and ground truth (right).
Figure 4. Top row: propose method with the shrinkage function (left) and
g [ j ],geo
g [ j ],sc
(right). Bottom row: SWT-MMSE (left) and Refined Lee Filter 9x9 (right).
Figure 5. Original Image (Top Left), proposed method with g [ j ],sc (top right), SWT-MMSE (bottom left) and Refined Lee 9x9 (bottom right). Images are displayed in db.
V. CONCLUSION In this paper, we proposed a new approach for the filtering of SAR images based on the stationary wavelet transform. The denoising of wavelet coefficients is achieved by a shrinkage function derived from a two component mixture model of the pdf of the normalized wavelet coefficients. This approach is particularly suitable for the filtering of correlated speckle and does not require any knowledge of the speckle second-order statistics. Different shrinkage functions are also proposed that are reinforcing consistent values and geometry across scales. Filtering results on artificial and real SAR images show a more suitable compromise between noise reduction and detail preservation as compared to standard filters. Noise reduction is stronger to what can be achieved by other wavelet based speckle reduction techniques that are assuming a white speckle noise. Among the four different shrinkage functions proposed, g [ j ],geo produces a speckle reduction similar to the Enhanced Lee filter (9x9) but with superior point target and edge preservation. The additional direction consistency rule (14) is not improving filtering results because point targets are altered. ACKNOWLEDGMENT This work has been supported in part by the NSERC of Canada (Discovery Grant) and the MDEIE of the “Gouvernement du Québec”. The authors wish to thank Elaine Rosenberg and Lisa Hollinger for their linguistic expertise. APPENDIX The following pseudo-code describes the estimation procedure for σn[ j ] and σe[ j ] (see section III.A). Initialization H ← histogram(M [ j ]I ) σˆn[ j ] ← arg max{H }/ 2 3 σ[ j ] ← ( ) M [ j ]I /( 2Γ(2)) 2 m ←1 α ← αmin While α < αmax σe[ j ](m ) ← (σ[ j ] − α(m )σˆn[ j ] )/(1 − α(m )) RMS (m ) ← computeRMS (H , p[ j ](r )) m ← m +1 α(m ) ← α(m − 1) + ∆α End σˆe[ j ] ← σe[ j ](arg min{RMS }) m
[j ]
πˆ
← α(arg min{RMS }) m
REFERENCES [1]
S. Foucher, G. B. Bénié, and J.-M. Boucher, Multiscale MAP Filtering of SAR images, IEEE Trans. on Image Proc., vol. 10, no. 1, pp. 49 –6, Jan. 2001. [2] A. Achim, P. Tsakalides and A. Bezerianos; “SAR Image Denoising via Bayesian Wavelet Shrinkage Based on Heavy-Tailed Modeling”, IEEE Trans. on Geosci. and Remote Sens., vol 41, no 8, pp. 1773 –1784, Aug. 2003. [3] Hua Xie, L. E. Pierce and F. T. Ulaby, “SAR Speckle Reduction using Wavelet Denoising and Markov Random Field Modeling”, IEEE Trans. on Geosci. and Remote Sens., vol. 40, no. 10, pp. 2196 –2212, Oct. 2002. [4] F. Argenti and L. Alparone, “Speckle Removal from SAR Images in the Undecimated Wavelet Domain”, IEEE Trans. on Geosci. and Remote Sens., vol. 40, no. 11, pp. 2363 –2212, Nov. 2002. [5] S. Foucher, J.-M. Boucher and G. B. Bénié. "Wavelet Filtering of Correlated Speckle" In AeroSense 2003 (SPIE #5108). Orlando, 21-25 April 2003. [6] J. Scharcanski; C. R. Jung, and R. T. Clarke, “Adaptive Image Denoising Using Scale and Space Consistency”, IEEE Trans. on Image Proc., vol. 11, no. 9, pp. 1092 –1101, Sept. 2002 [7] A. Papoulis, Probability, “Random Variables, and Stochastic Processes,” McGraw-Hill, 1991. [8] D. Blacknell, “New method for the simulation of correlated k-distributed clutter,” IEE Proceedings, radar, Sonar Navigation 141, p. 45, Feb. 1994. [9] J.-S. Lee, M.R. Grunes, and G. de Grandi, “Polarimetric SAR Speckle Filtering and Its Implication for Classification,” IEEE Trans. on Geosci. and Remote Sens., vol. 37, no. 5, pp. 2363-2373, 1999. [10] S. Foucher, G. Farage, and G. B. Bénié, “SAR Image Filtering Based on the Stationary Contourlet Transform”, IGARSS, Denver, pp. 4021-4024 Jul. 2006. [11] G. Farage, S. Foucher, and G. B. Bénié, “Comparison of PolSAR Speckle Filtering Techniques”, IGARSS, Denver, pp. 1760-1763, Jul. 2006.
