Edge-Preserving Unscented Kalman Filter for ... - Semantic Scholar

Report 2 Downloads 152 Views
Edge-Preserving Unscented Kalman Filter for Speckle Reduction G.R.K.S. Subrahmanyam1 A.N. Rajagopalan1 R. Aravind 1 Gerhard Rigoll2 1 Department of Electrical Engineering Indian Institute of Technology Madras, Chennai 600 036, India 2 Lehrstuhl f¨ ur Mensch-Maschine-Kommunikation Technische Universit¨ at M¨unchen, 80333 M¨unchen, Germany {ee03d009,raju,aravind}@ee.iitm.ac.in, [email protected]

Abstract We propose a recursive spatial-domain speckle reduction algorithm for synthetic aperture radar (SAR) imagery based on the unscented Kalman filter (UKF) with a discontinuity-adaptive Markov random field (DAMRF) prior. The capability of the UKF in handling speckle noise and the feature preservation ability of the DAMRF model are explored within a unified framework through importance sampling.

1 Introduction Speckle noise severely impedes automatic scene segmentation and interpretation, and limits the resolution of SAR images as well as their utility. If s represents the original image and v is speckle noise, then the degraded observation y is given by the relation y(m, n) = s(m, n) · v(m, n)

(1)

Noise v is assumed to be independent of s with unit mean and variance σv2 . The multiplicative nature of speckle complicates the noise filtering process in SAR. A speckle suppression filter is expected to effectively filter homogeneous areas, retain image texture and edges, and preserve features (both linear and pointtype). In the Lee filter [1], the multiplicative model is approximated by a linear combination of the local mean and the observed pixel, and a minimum mean square error estimator is applied to determine the weighting constant. The Frost filter [2] is an adaptive and exponentially-weighted averaging filter based on the coefficient of variation c which is the ratio of the local standard deviation to the local mean of the degraded image. The enhanced Lee and Frost filters proposed by Lopes et al. [3] divide an image into homogeneous,

978-1-4244-2175-6/08/$25.00 ©2008 IEEE

heterogeneous areas, and isolated point targets based on the coefficient of variation c (low, intermediate, and high, respectively). In [4], speckle is suppressed by employing an adaptive block-Kalman filter, which relies on block-varying autoregressive (AR) parameters of the original image. An adaptive maximum a posteriori (MAP) estimator with a heavy-tailed Rayleigh model has been suggested in [5]. A wavelet despeckling method based on Bayesian shrinkage has been proposed in [6]. Argenti et. al. [7] propose despeckling in the undecimated wavelet domain using a space-varying generalized Gaussian distribution for the wavelet coefficients. Adaptive filters [1, 2, 3] are simple but tend to oversmooth image texture and suffer from ineffective denoising around edges. The adaptive block-Kalman filter [4] requires the AR coefficients of the noise-free image to be known accurately. Recently proposed wavelet domain [6] and MAP approaches [5] have superior performance but require explicit optimization and/or parameter estimation. We propose a novel method based on the unscented Kalman filter (UKF) [8, 9] to suppress speckle noise in SAR images. Multiplicative speckle noise is handled with the measurement equation of the UKF. A small set of sigma points which can capture and propagate the first two statistics through the measurement model are employed in a recursive framework for image estimation. A non-Gaussian discontinuity adaptive Markov random field (DAMRF) prior is incorporated into the UKF framework through a Monte Carlo approach known as importance sampling (IS). Hence, we refer to our method as ISUKF.

2

Principle of UKF

The UKF is based on the unscented transformation (UT) which is a deterministic sampling approach for calculating the statistics of a random variable which

undergoes a transformation. Consider propagating an nx -dimensional random variable x through a function g : Rnx →Rny to generate y = g(x). Assume x has mean x and covariance Px . To calculate the statistics (first two moments) of y using the scaled unscented transformation (SUT), we proceed as follows: First, a set of 2nx + 1 weighted samples or sigma points are deterministically chosen so that they exactly capture the true mean and covariance of the prior random variable x. A selection scheme that satisfies this requirement [9] is p X0 = x; Xi = x + ( (nx + λ)Px )i , i = 1, ...., nx ; p Xi = x − ( (nx + λ)Px )i , i = nx + 1, ...., 2nx ; λ (m) λ = α2UT (nx + κ) − nx w0 = (nx + λ) (c)

(m)

+ (1 − α2UT + βUT ); 1 , i = 1, ...., 2nx ; (2) = 2(nx + λ)

w0 = w0 (m)

wi

(c)

= wi

Here, scaling parameters αUT p, βUT and κ are chosen as suggested in [9], while ( (nxp + λ)Px )i is the ith column of the matrix square root ( (nx + λ)Px )i [8]. Superscripts m or c on the weight wi refer to their use in mean or covariance calculation, respectively. The sigma point matrix X is formed by stacking the (2nx + 1) sigma points Xi (column wise) and hence has dimension nx × (2nx + 1). Simple matrix calculations confirm that the first two moments of these sigma points are exactly the same as that of the prior [8]. Each sigma point is passed through the nonlinearity to get Yi = g(Xi ), i = 0, 1, ....., 2nx. The mean and covariance of y are estimated as y=

2nx X i=0

(m)

wi

Yi ,

Py =

2nx X

   1 η 2 (s(m, n), sb) exp −γ log 1 + Z γ

(4)

= where P energy term η 2 (s(m, n), sb) 1 (s(m, n) − s b (m − i, n − j))2 and (i,j)∈R 2ρ2m,n (i, j) ∈ R the neighborhood of (m, n). Since it is analytically not possible to estimate the mean and variance of this non-Gaussian pdf, we resort to a Monte Carlo method known as importance sampling (IS) [12] which is a method for determining the estimates under a (difficult to sample) target pdf, say p(z), provided its functional form is known up to a multiplication constant [12]. Consider a distribution q(z), a reasonable approximation to p(z), which is also known up to a multiplicative constant, is easy to sample, and is such that the (non-zero) support of q(z) includes the support of p(z). Such a density q(z) is called the sampler density (referred to as importance function). This is based on the Monte Carlo approximation   Z Z p(z) Ep [g(z)] = g(z)p(z)dz = g(z) q(z)dz = q(z)      L 1X p(z) p(z (l) ) ≈ (5) Eq g(z) g(z (l) ) q(z) L i=1 q(z (l) ) where L samples {z (l) } are drawn from the distribution given by q(z). The selection of the sampler density q is crucial in determining the accuracy of the IS estimates.

4

Speckle Reduction using ISUKF

In this section, we present the steps in the proposed approach for estimating each pixel sequentially.

(c)

wi (Yi −y)(Yi −y)T .

i=0

(3) These estimates of the moments of y are accurate to the second-order (and third-order for symmetric priors) of the Taylor series expansion for any nonlinear function g(x), since the prior sigma points completely capture the prior distribution upto the second moment [8].

3 Moments of Non-Gaussian Prior The well-known Gaussian MRF model is effective in modeling smooth data, but it penalizes edges. Geman and Geman [10] proposed the use of line fields model for preserving edges but at the cost of rendering the potential function non-differentiable. Following Li [11], we model image s by a continuous discontinuity adaptive MRF (DAMRF) with conditional pdf P (s(m, n)|b s(m − i, n − j)) =

1. At pixel (m, n), using the DAMRF model and the past pixel estimates, we construct the state conditional pdf as p (s(m, n)|b s(m − i, n − j)) =    η 2 (s(m, n),bs) exp −γ log 1 + γ

(6)

Here, P we define η 2 (s(m, n),bs) = 1 2 (s(m, n) − s b (m − i, n − j)) 2 (i,j) 2ρm,n and for a first-order NSHP neighborhood (i, j) ∈ {(0, 1), (1, 0), (1, 1), (1, −1)}. For effective smoothing of speckle and preservation of fine features, we vary the scale parameter as ρ2m,n = ζµ1.5 s , where ζ is user-specified and µs is the mean of the first-order past NSHP pixels. 2. We employ importance sampling (IS) technique to estimate the first-two moments of the above

conditional pdf as follows. We draw samples {z l }, l = 1, 2, ..., L from a Cauchy sampler as it has a heavy tail. To facilitate effective IS, the location µq of the Cauchy sampler density q(z) is varied as µq = µs and the scale is chosen sufficiently high to include the support of p in q. The samples are p(z l ) unbiased using the correction weights wl = q(z l ) (following Eq. 5). The mean µ bp and variance σ bp2 of p are computed as PL PL (l) (l) −µ bp )2 2 l=1 wl z l=1 wl (z µ bp = P σ b = . P p L L l=1 wl l=1 wl 3. The above predicted state statistics are augmented with measurement noise statistics and the sigma point selection scheme (Eq. 2) is applied: xa(m,n)/(m,n−1) = [b µp 1]T   2 σ bp 0 a P(m,n)/(m,n−1) = 0 σv2 h Xa(m,n)/(m,n−1) = xa(m,n)/(m,n−1) i q xa(m,n)/(m,n−1) ± (na + λ)Pa(m,n)/(m,n−1)

Here, Xa is the augmented sigma point matrix (since na = 2, following section (2), Xa is of dimension 2 × 5) and Xa = [(Xx ) (Xv )]T with Xx and Xv sigma points corresponding to state and measurement noise, respectively, each being 1 × 5. Now, we independently propagate each of these sigma points through the measurement model to predict the sigma points corresponding to measurement as: Y(m,n)/(m,n−1) = Xx . ∗ Xv Using sigma points Xx and Y, we estimate y (m,n)/(m,n−1) , Pxy and Pyy by Eq. (3).

(a)

(b)

(c)

(d)

Figure 1. (a) Original. (b) Noisy. Output of (c) Frost (F OM = 0.94, S/M SE = 20.51 dB), and (d) the proposed method (F OM = 0.98, S/M SE = 22.83 dB).

Consider the edge image shown in Fig. 1(a). This image is degraded with simulated 2-Look Gammadistributed speckle noise (Fig. 1(b)). The outputs of the Frost [2] and the proposed ISUKF method are shown in Figs. 1(c), and (d), respectively. The image in Fig. 1 (c) has a grainy appearance in the white uniform region due to residual speckle. The proposed method is not only effective in despeckling but also well preserves the transition between the dark and bright regions. This is also reflected in the S/M SE and F OM values.

The Kalman gain K(m,n) = Pxy P−1 yy . We update the predicted mean as: sb(m, n) = s(m,n)|(m,n−1) + K(m,n) (y(m,n) − y (m,n)/(m,n−1) )

5 Experimental Results For ISUKF, we set αUT = 1, βUT = 0, κ = 0, ζ = 0.01, and L = 100. Quantitative evaluation is done with Signal to Mean Square Error ratio (S/M SE), Equivalent Number of Looks (EN L) which measures the speckle suppression in a homogeneous area of the image, and Figure of Merit (0 ≤ F OM ≤ 1) which evaluates edge rendition. Higher the values of these metrics, better the performance.

(a)

(b)

Figure 2. (a) The Horse track image. (b) Output of our method.

Next, we considered the case of a 1-look ’Horse track’ image that is affected by real speckle noise (Fig. 2(a)). The image estimated by the proposed ISUKF

(a)

(b)

(c)

Figure 3. (a) Bedfordshire image in Southeast England (EN L = 12.16). Result using (b) edgebased wavelet filter [6] (EN L = 42.87), and (c) the ISUKF method (EN L = 46.91).

method is shown in Fig. 2(b). Our method is quite effective in smoothing speckle over uniform regions (such as the top-left and bottom-right gray regions); yet the edges are sharp and clear. Figure 3(a) shows a 2-look SAR image with real speckle. The image estimated using a recently proposed wavelet filter is reproduced from [6] in Fig. 3(b). The output of the ISUKF filter (Fig. 3(c)) is comparatively superior. We note that the edges are recovered without any blurring effects while the homogeneous regions are almost free from noise. The isolated point targets (such as the two white spots in the dark homogeneous region on the left) become strikingly visible in Fig. 3(c).

6 Conclusions We proposed a novel and effective UKF-based speckle reduction algorithm. A small set of sigma points was used to propagate the desired statistics through the UKF filter equations. A non-Gaussian prior was incorporated within the recursive UKF framework through importance sampling for feature preservation. On a Pentium-4 PC with 256 MB RAM and for a 200 × 200 image, our method takes less than 20 seconds to execute. Acknowledgment: The second author is grateful to the Humboldt Foundation for its support.

References [1] J. S. Lee, “Digital image enhancement and noise filtering by use of local statistics”, IEEE Trans. Pattern Anal. Machine Intell., vol. 2, pp. 165-168, 1980. [2] V.S. Frost, J.A. Stiles, K.S. Shanmugan, and J.C. Holtzman, “A model for Radar images and its application to

adaptive digital filtering of multiplicative noise”, IEEE Trans. Pattern Anal. Machine Intell., vol. 4, pp. 157-166, 1982. [3] A. Lopes, R. Touzi, and E. Nezzy, “Adaptive speckle filters and scene heterogeneity,” IEEE Trans. Geoscience and Remote Sensing, vol. 28, pp. 992-1000, 1990. [4] M. R. Azimi-Sadjadi and S. Bannour, “2D adaptive block Kalman filtering of SAR imagery”, IEEE Trans. Geoscience and Remote Sensing, vol. 29. pp. 742-753, 1991. [5] A. Achim, E. E. Kuruoglu and J. Zerubia, “SAR image filtering based on the heavy-tailed Rayleigh model”, IEEE Trans. Image Processing, vol. 15, pp. 2686-2693, 2006. [6] M. Dai, C. Peng, A. K. Chan, and D. Loguinov, “Bayesian wavelet shrinkage with edge detection for SAR image despeckling”, IEEE Trans. Geoscience and Remote Sensing, vol. 42, pp. 1642-1648, 2004. [7] F. Argenti, T. Bianchi and L. Alparone, “Multiresolution MAP despeckling of SAR images based on locally adaptive generalized Gaussian pdf modeling”, IEEE Trans. Image Processing, vol. 15, pp. 3385-3399, 2006. [8] S. J. Julier and J. K. Uhlmann. “A general method for approximating nonlinear transformations of probability distributions”. Technical report, RRG, Dept. of Engineering Science, University of Oxford, Nov 1996. [9] S. J. Julier and J. K. Uhlmann, “A new extension of the Kalman filter to nonlinear systems”, Proc. of AeroSense: The 11th Intnl. Symposium on Aerospace/Defense Sensing, Simulation and Controls, Orlando, Florida, 1997. [10] S. Geman and D. Geman, “Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images”, IEEE Trans. Pattern Anal. Machine Intell., vol. 6, pp. 721741, 1984. [11] S. Z. Li, Markov Random Field Modeling in Computer Vision, Springer, 1995. [12] D. J. C. Mackay, “Introduction to Monte Carlo methods”, In M. I. Jordan, Ed., “Learning in graphical models”, NATO Science Series, pp. 175–204, Kluwer, 1998.