Feature Article
Importance Sampling-Based Unscented Kalman Filter for Film-Grain Noise Removal G.R.K. Sai Subrahmanyam, A.N. Rajagopalan, and Rangarajan Aravind Indian Institute of Technology Madras Photographic film contains film-grain noise that translates to multiplicative, non-Gaussian noise in the exposure domain. A method based on the unscented Kalman filter can suppress this noise while simultaneously preserving edge information.
52
P
hotographic film is a widely used image-recording medium. Negatives contain tiny crystals of silver halide salts, which are light sensitive. When such a film is developed, these crystals turn into tiny filaments of metallic silver conventionally called grains. These randomly patterned grains, which constitute filmgrain noise, are often visibly objectionable in photographic prints. While film grain can enhance the feel of a motion picture,1 it renders the task of compression (transmission) difficult because of the noise. For the multimedia industry, reducing film-grain noise is important for scanning, digitization, and efficient storage.2 There is a well-known nonlinear relationship between the incoming light intensity (exposure) and the silver density deposited on the film. The domain in which the actual film formation takes place is referred to as the density domain and the domain in which the developed image is available for visual consumption is referred to as the exposure domain.
In the density domain, an image affected by film-grain noise can be modeled as a logarithmic nonlinearity of the original scene corrupted by additive white Gaussian noise.3,4 That is, rd(m,n) 5 alog10(s(m,n)) + b + v(m,n) where s(m,n) is the original scene, rd(m,n) is the degraded observation in the density domain, and v(m, n) denotes additive white Gaussian noise with zero mean and variance s2v . Parameters a and b are the slope and offset derived from the film’s D 2 logE curve. Alternatively, in the exposure domain, we can write re ðm,nÞ ~ sðm,nÞ 10vðm,nÞ=a
ð1Þ
but the noise becomes multiplicative and non-Gaussian. In Equation 1, re(m,n) is the degraded exposure domain image and is related to the density domain observation image as rd(m,n) 5 alog10(re(m,n)) + b. The unscented Kalman filter (UKF) can incorporate multiplicative, non-Gaussian noise into the estimation procedure. Demonstrating this capability in the 2D domain, we propose an autoregressive UKF for film-grain denoising. When modeled as an AR process, the image prior has limited capability to handle edges. We propose a novel methodology that reduces film-grain noise but preserves the edge information by judiciously incorporating a non-Gaussian prior within the (UKF) framework through importance sampling. We achieve this by formulating a discontinuityadaptive Markov random field (DAMRF) prior suited for recursive prediction and by employing importance sampling with a space-varying and heavy-tailed proposal density to estimate the DAMRF statistics.
State-space formulation and unscented transformation The UKF framework is similar to that of the Kalman filter. Like the Kalman filter, the UKF is based on state-space formulation but differs in the way it represents and propagates Gaussianity through system dynamics. If the image is modeled as a 2D autoregressive (AR) process, the corresponding AR equation written in state transition form is xðm,nÞ ~ f xðm,n{1Þ ,uðm,nÞ ~ Fxðm,n{1Þ z uðm,nÞ
1070-986X/08/$25.00 G 2008 IEEE Published by the IEEE Computer Society Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
Here, x(m,n) is the state vector with dimension nx. If we consider a first-order, causal, nonsymmetric half plane (NSHP) support, then nx 5 3, x(m,n) 5 [s(m,n), s(m, n 2 1), s(m 2 1,n)]T and x(m,n21) 5 [s(m,n 2 1),s(m 2 1, n),s(m 2 1, n 2 1)]T where (m, n) denotes the current pixel position, and superscript T refers to transpose. The state transition matrix F contains the AR coefficients of the image. Vector u(m,n) 5 (u(m,n) 0 0)T where u(m,n) is state noise with variance s2u and dimension nu(5 1). For the film-grain noise problem, the measurement equation in exposure domain (see Equation 1) is of the form
yðm,nÞ ~ h xðm,nÞ ,vðm,nÞ ~ Hxðm,nÞ 10vðm,nÞ =a
i
z 1{ z bUT pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi { ðnx z lÞP , i ~ nx z1, . . . ,2nx ; Xi ~ x
ðc Þ w0
~
ðmÞ w0
~
ðcÞ wi
a2UT
i
ðmÞ wi
1 ~ , i ~ 1, . . . ,2nx ð2Þ 2ðnx z lÞ
Here, l ~ a2UT ðnx z kÞ { nx . Parameter aUT controls the spread of the sigma points around x ¯ and is usually set to a value between 0 and 1. The term bUT is used to incorporate prior knowledge about x. Note pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðnx z lÞP i represents the ith row (or that column) of the nx 3 nx-matrix square root of (nx + l)P.7 The addition and subtraction of these vectors from the mean value of the prior results in 2nx symmetric sigma points. Along with the prior mean X0 we have 2nx + 1 sigma points. Subscript i on X refers to the ith sigma point. The superscripts m and c on the weights refer to their use in mean and covariance calculations, respectively. The ðmÞ
weight wi associated with the ith point P 2nx ðmÞ satisfies ~ 1. Simple matrix calcui~0 wi lations can confirm that the first two moments of these sigma points are exactly the same as that of the prior. We instantiate each sigma point through the true nonlinearity to yield Yi 5 g(Xi),i 5 0,1,…,2nx. We estimate the mean and covariance of y as
y~
X 2n
ðmÞ x i~0 wi Yi
and Py ~
X 2n
ðcÞ x i~0 wi ðYi
{ yÞðYi { yÞT
Positive semidefiniteness of the predicted covariance matrix Py is guaranteed by choosing k $ 0.6,7 These estimates of the moments of y are accurate to the second-order (thirdorder for symmetric priors) of the Taylor series expansion for any nonlinear function g(x), because the prior sigma points completely capture the prior distribution up to the second moment.7
April–June 2008
where y(m,n) 5 re(m,n) is a scalar measurement (with dimension ny 5 1), v(m, n) 5 v(m, n) is measurement noise of dimension nv(5 1), and H 5 [1 0 0]. For the film-grain noise problem, the measurement model h is nonlinear. Hence, we resort to a nonlinear counterpart of the Kalman filter, namely the UKF. The UKF is a straightforward application of the unscented transformation (UT) for recursive minimum mean-square-error estimation.5,6 The UT is founded on the notion that with a fixed number of parameters, it’s easier to approximate a Gaussian distribution than to approximate an arbitrary nonlinear function or transformation.5,7 The UT is a deterministic sampling approach for calculating the statistics of a random variable that undergoes a nonlinear transformation. Consider propagating the random variable x through a nonlinear function g : Rnx ?Rny to generate y 5 g(x). In our problem, the function g is equivalent to f for the state transformation, and h for the measurement transformation. Assume x has mean x ¯ and covariance P. To calculate the first two moments of y using the scaled UT, we proceed as follows: First, we deterministically choose a set of 2nx + 1 weighted samples or sigma points so that they exactly capture the true mean and covariance of the prior random variable x. A selection scheme that satisfies this requirement,5,6 is as follows:
l ðnx z lÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi z ðnx z lÞP , i ~ 1, . . . ,nx ; Xi ~ x ðmÞ
; w0 ~ X0 ~ x
53 Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
Auto-regressive UKF
where Xxðm,n{1Þ
In this section, we propose a recursive filter, based on the UKF, that accounts for the image sensor nonlinearity in the observation model. The UKF, based on the UT, provides a mathematically tractable way to propagate the first two moments even in the presence of multiplicative and non-Gaussian noise. We achieve this by augmenting the state-random variable with the (state and observation) noise random variables and applying UT in the prediction of state and observation.5 We apply the scaled UT sigma point selection scheme to the new augmented vector h iT xaðm,nÞ ~ xTðm,nÞ uðm,nÞ vðm,nÞ with dimen-
formed from the first nx rows of Xaðm,n{1Þ ,
sion na to calculate the corresponding augmented sigma point matrix Xaðm,nÞ of dimension na 3 (2na + 1). Here na 5 nx + nu + nv. (For example, for a first-order NSHP support nx ~ 3,X aðm,nÞ will be of size 5 3 11 and each sigma point is of dimension 5). We now propose the auto-regressive UKF (ARUKF), which sequentially estimates the intensity of the original image at each pixel (m, n) in a raster scan order from left to right. The filter updates the mean and covariance of the Gaussian approximation to the posterior distribution of the state as we describe in the following steps. We start by initializing state statistics x ¯0 5 E[x0];P0 = E[(x0 2 x ¯0)(x0 2 x ¯0)T]. Step 1 augments state mean and covariance with noise statistics: h
iT
For step 2, we calculate the augmented sigma point matrix (see Equation 2):
IEEE MultiMedia
Xaðm,n{1Þ ~ h qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i aðm,n{1Þ + ðna z lÞPaðm,n{1Þ aðm,n{1Þ x x In step 3, the state sigma points and state noise sigma points are propagated through the AR state model for prediction: Xxðm,nÞ=ðm,n{1Þ ~ f Xxðm,n{1Þ ,Xuðm,n{1Þ ~ F Xxðm,n{1Þ
while Xuðm,n{1Þ is formed from the nu rows that immediately follow the first nx rows in Xaðm,n{1Þ . For step 4, we use the predicted state sigma points to determine the predicted mean and covariances: ðm,nÞ=ðm,n{1Þ ~ x
X 2n
ðmÞ x a i~0 wi Xi,ðm,nÞ=ðm,n{1Þ
X 2n ðcÞ a Pðm,nÞ=ðm,n{1Þ ~ i~0 wi h i ðm,nÞ=ðm,n{1Þ | Xxi,ðm,nÞ=ðm,n{1Þ { x h iT ðm,nÞ=ðm,n{1Þ | Xxi,ðm,nÞ=ðm,n{1Þ { x
where i represents the ith sigma point that is the ith column of Xxðm,nÞ=ðm,n{1Þ . In step 5, the sigma points corresponding to the measurement are predicted by propagating the (predicted) state sigma points and measurement noise sigma points through the filmgrain nonlinearity in the exposure domain: Yðm,nÞ=ðm,n{1Þ ~ h Xxðm,nÞ=ðm,n{1Þ ,Xvðm,n{1Þ v ~ HXxðm,nÞ=ðm,n{1Þ: 10: Xðm,n{1Þ a Here, Xvðm,n{1Þ is formed from the last nv rows of Xaðm,n{1Þ . Operations .* and .() represent element-by-element operations (multiplication and raised power, respectively). This results in measurement sigma point matrix Y(m, n)/(m, n21) of dimension ny 3 (2na + 1). In step 6, we estimate the statistics required for updating using the predicted sigma points:
aðm,n{1Þ ~ x Tðm,n{1Þ 0 0 x 2 3 Pðm,n{1Þ 0 0 6 7 Paðm,n{1Þ ~ 6 0 s2u 0 7 4 5 2 0 0 sv
h z XuT ðm,n{1Þ
consists of sigma points
0T
0T
iT
X 2n ðmÞ a yðm,nÞ=ðm,n{1Þ ~ i~0 wi Yi,ðm,nÞ=ðm,n{1Þ i X 2n ðcÞ h a Pyy ~ Yi,ðm,nÞ=ðm,n{1Þ { yðm,nÞ=ðm,n{1Þ i~0 wi h iT yðm,nÞ=ðm,n{1Þ | Yi,ðm,nÞ=ðm,n{1Þ { i X 2n ðcÞ h x a Pxy ~ w X { x ð m,n Þ= ð m,n{1 Þ i,ðm,nÞ=ðm,n{1Þ i~0 i h iT yðm,nÞ=ðm,n{1Þ | Yi,ðm,nÞ=ðm,n{1Þ { where Yi, (m, n)/(m, n21) represents the ith sigma point, which is the ith column of Y(m, n)/(m, n21).
54 Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
Step 7 performs updates, as in the Kalman filter: The Kalman gain Kðm,nÞ ~ Pxy P{1 yy ðm,nÞ=ðm,n{1Þ ðm,nÞ ~ x x z Kðm,nÞ yðm,nÞ { yðm,nÞ=ðm,n{1Þ Pðm,nÞ ~ Pðm,nÞ=ðm,n{1Þ { Kðm,nÞ Pyy KTðm,nÞ Estimated pixel intensity at (m, n) is ˆs(m,n) 5x ¯(m,n)(1), the first component of x ¯(m,n).
Non-Gaussian MRF prior
Pðsðm,nÞ=^sðm { i,n { jÞÞ ~ 1 g2 ðsðm,nÞ,^sÞ exp {clog 1 z Z c
ð3Þ
where Z is a normalization constant and g2 ðsðm,nÞ,^sÞ ~ 1
X ðsðm,nÞ { ^sðm { i,n { jÞÞ2 ði2 z j2 Þ
r2ðm,nÞ ði,jÞ[R
Here, (i, j) M R 5 {(i, j)I(1 # i # M1, 2M1 # j # M1) < (0, 1 # i # M1)}, and M1 is the order of the NSHP support. Figure 1a (next page) illustrates NSHP support. Parameter r2ðm,nÞ controls the allowable variation of the current pixel with respect to its neighboring values.15 Here, the potential function that accounts for the irregularities in the solution is given by 2 c log 1 z gc . The plot in Figure 1b shows the degree of penalty dictated by the potential function. For this MRF model, the smoothing strength of the regularizer increases monotonically as g increases within the convex band pffiffiffi pffiffiffi Bc ~ { c, c . Outside the band, smoothing decreases and becomes zero as g R ‘. Because this feature preserves image discontinuities, it’s called a DAMRF model.13
Importance sampling To incorporate the non-Gaussian prior given by the DAMRF model (see Equation 3) into the UKF framework, we need to estimate the mean and covariance of this prior. Because it is analytically impossible to compute these moments, we resort to a Monte Carlo method known as importance sampling.16 This technique determines a target PDF’s estimates using the samples from an easy-to-sample approximation. Consider a PDF p(z) that is known up to a multiplicative constant but is difficult to esti-
April–June 2008
A homogeneous and linear AR model imposes a strong constraint on state evolution, which can hinder true transitions at image edges. Other research has addressed this issue. For example, the research of Woods et al. uses spatially varying 2D AR parameters, estimated by windowing the observed image.8 Jeng and Woods propose inhomogeneous AR image models using local statistics.9 Kadaba et al. observe that Bernouli-Laplacian and Cauchy distribution are a better fit than the Gaussian for error residuals.10 Jeng and Woods use a compound Gauss MRF to model images.11 Ceccarelli proposed edge-preserving MRFs.12 We base our approach to incorporating a non-Gaussian prior within our recursive estimation scheme on the fact that statistical dependence incorporated using an MRF model provides more flexibility in incorporating contextual constraints through a suitably derived prior. An MRF possesses a Markovian property—that is, the value of a pixel depends only on the values of its neighboring pixels.13 It’s important to realize that at discontinuities or edges, the notion of neighborhood dependency should be relaxed. Incorporating the smoothness constraint while simultaneously preserving the edges is a challenging task. For simplicity and analytical tractability, it is normal practice to model the original image s by a homogeneous Gaussian MRF. The issue with a homogeneous Gaussian MRF is that it imposes smoothness constraints everywhere, which inevitably leads to oversmoothing edges. A better way to handle the situation is to use a non-Gaussian conditional probability density function (PDF) to model the original image. Geman and Geman propose the use of line fields in which the smoothness constraint is switched off at points where the magnitude
of the signal derivative exceeds a certain threshold.14 However, the resulting potential function becomes discontinuous. Li proposes a continuous adaptive regularizer model in which, unlike the line-field model, the interaction decreases as the derivative magnitude becomes larger and is completely turned off at infinite magnitude of the signal derivative.13 Following Li, we model s as a non-Gaussian MRF with conditional PDF given by
55 Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
p ðz Þ qðzÞdz q ðz Þ !
L pzðlÞ pðzÞ 1X ðlÞ ~ E q g ðz Þ g z & qðzÞ L l~1 qðzðlÞ Þ
Figure 1. (a)
E p ½ g ðz Þ ~
Nonsymmetric halfplane support at a pixel (m, n). (b) The degree of penalty with distance imposed by
ð
g ðzÞpðzÞdz ~
ð
g ðzÞ
Here, z(l) represents samples drawn from q(z). The basic idea behind importance sampling is to choose a distribution that emphasizes certain values of the input random variables that have more impact on the parameter being estimated than others by sampling them more frequently. Our aim is to determine the estimates under the PDF p, using samples from q. Because it’s difficult to draw samples from L the target PDF p, we draw L samples, zðlÞ l~1 from the sampler PDF q—that is, instead of choosing samples from a uniform distribution, they are now chosen from a distribution q that concentrates on those points where the function p is large. When we use samples from q to determine any estimates under p, in the regions where q is greater than p, the estimates are over-represented. In the regions where q is less than p, they are underrepresented. To account for this, we use correction weights
the discontinuityadaptive Markov random field model. (c) Importance sampling: p is the target probability density function whose moments are to be estimated, while q is the (Cauchy) sampler density.
p zðlÞ w ~ qðzðlÞ Þ l
in determining the estimates under p. For example, to find the mean of the distribution p we use
IEEE MultiMedia
^ mp ~
mate its moments. From the functional form, we can estimate its support. 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 the (non-zero) support of q(z) includes the support of p(z). Such a density q(z) is called the sampler density. Figure 1c shows a typical plot that illustrates the target PDF p and sampler distribution q. Mathematically, we can explain the principle of importance sampling as follows. Suppose the density q(z) roughly approximates the density of interest p(z), then the expected value of a function g of a random variable z is
PL wl zðlÞ Pl~1 L wl l~1
As L R ‘, the estimate mˆp tends to the actual mean value of p. In our problem, the DAMRF distribution corresponds to p. We choose the Cauchy distribution for the importance function q, as it yields better estimates (due to its heavy tailed nature) compared to the Gaussian sampler.16
Edge-preserving UKF We now extend the basic structure of the UKF to include a non-Gaussian prior. In the proposed strategy, we use knowledge of the conditional PDF to capture the image’s local physical characteristics. Doing so implicitly generalizes the state transition equation and
56 Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
doesn’t restrict it to being linear. We base our proposed extension on the observation that the Kalman filter variants differ only in state vector propagation but have same update equations.5 Specifically, we formulate the state conditional density of the pixel to be estimated (given its neighbors) based on the DAMRF model discussed previously. We estimate the mean and variance of the non-Gaussian conditional prior, which constitutes the predicted statistics of the state, by importance sampling. On the basis of these and known measurement noise statistics, we determine predicted sigma points and propagate them through the nonlinearity. The approach takes the sigma points to the update stage of the UKF and determines the final image estimate. The steps involved in the proposed importance sampling UKF (ISUKF) method follows. For step 1, at each pixel, construct the state conditional PDF using the past pixels in the NSHP support, and the values of r2ðm,nÞ and c in the DAMRF model. That is, Pðsðm,nÞ=^sðm { i,m { jÞÞ g2 ðsðm,nÞ,^sÞ ~ exp {c log 1 z c If M1 5 1, then the neighbors considered in the first-order NSHP support are ˆs(m,n 2 1), ˆs(m 2 1,n), ˆs(m 2 1,n 2 1), and ˆs(m 2 1,n + 1), which already are estimated. Because the parameter r2ðm,nÞ of the DAMRF model charac-
p zðlÞ l w ~ qðzðlÞ Þ
^p ~ m
PL wl zðlÞ Pl~1 L wl l~1
and
^2p s
PL ~
l~1 w
l
^p zðlÞ { m
PL
l~1 w
2
l
In step 3, we directly use these estimates to predict the one-step-ahead sigma points. These estimates for the mean and covariance are augmented to deal with the nonadditive noise, as the following shows: ðm,nÞ=ðm,n{1Þ ~ ^ mp , x ^2p , Pðm,nÞ=ðm,n{1Þ ~ s h Tðm,nÞ=ðm,n{1Þ aðm,nÞ=ðm,n{1Þ ~ x x
Paðm,nÞ=ðm,n{1Þ
iT 0
h iT ~ ^ mp 0 , and " # Pðm,nÞ=ðm,n{1Þ 0 ~ 0 s2v 2 3 ^2p 0 s 5 ~4 0 s2v
Because we base prediction directly on the prior PDF, we need not augment with the state noise. We use these augmented mean and covariances to determine the sigma point matrix (see Equation 2): Xaðm,nÞ=ðm,n { 1Þ ~ h aðm,nÞ=ðm,n { 1Þ x qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i aðm,nÞ=ðm,n { 1Þ + ðna z lÞPaðm,nÞ=ðm,n { 1Þ x
Then we propagate each of the 2na + 1 sigma points corresponding to state and measurement noise through the film-grain nonlinearity. Yðm,nÞ=m,n{1Þ ~ h Xxðm,nÞ=m,n{1Þ ,Xvðm,n{1Þ : 10: Xvðm,n{1Þ a ~ Xx ðm,nÞ=m,n{1Þ
April–June 2008
terizes the local dependency, we set it to P(m,n21), the covariance of the previous state. For step 2, obtain the mean and covariance of the previous PDF using importance sam L pling. Draw samples zðlÞ l~1 from a Cauchy sampler. Cauchy distributed samples with location k1 and scale s1 can be generated using k1 + s1 tan(p(u 2 0.5)), where u is uniformly distributed from 0 to 1. To facilitate a close match with the DAMRF PDF p at each pixel, the location of the sampler density q(z) is varied as mq ~ 14 ð^sðm,n { 1Þ z ^sðm { 1,nÞz ^sðm { 1,n { 1Þ z ^sðm { 1,n z 1ÞÞ and the scale is chosen so that it is large enough to include the support of p in q. The samples are weighted by correction weights
^2p of p are The mean mˆp and variance s computed as
57 Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
Figure 2. (a) Image of a
where X xðm,nÞ=m,n{1Þ and X vðm,n{1Þ are formed X aðm,nÞ=m,n{1Þ
house. (b) Degraded
from the first nx rows of
image (s2v ~ 0:05). (c)
rows following nx, respectively. Here .* and.() are element-by-element operations. Finally, we estimate the statistics of y (from the sigma points) and update following steps 6 and 7 as in the auto-regressive UKF. We derive the estimated pixel intensity by ˆs(m,n) 5 x ¯(m,n). Thus, ˆs yields the estimated image.
Result using a modified Wiener filter (improvement in signal-to-noise ratio 5 2.51 dB), (d) particle filter (ISNR 5 3.48 dB), (e) autoregressive unscented Kalman filter (ISNR 5 3.40 dB), and (f) importance sampling unscented Kalman filter (ISNR 5 4.55 dB).
and the nv
Experimental results Here, we present results obtained using the proposed ARUKF and ISUKF methods and compare their performance with the standard modified Wiener filter (MWF)4 and with the recently proposed particle filter (PF) method.17 The c and r(m,n) parameters influence the DAMRF model. We have found experimentally that c 5 1.8 yields the best performance. From the D 2 logE characteristics of the film, we
chose a 5 5 for all experiments. We set the parameters as aUT 5 1, bUT 5 0, and k 5 1. The UKF is robust to variations in these parameters. Simulated degradation For a quantitative comparison, we use the improvement in signal-to-noise ratio (ISNR), defined as P ISNR ~ 10log10 P
m,n ðre ðm,nÞ
{ sðm,nÞÞ2
sðm,nÞ m,n ð^
{ sðm,nÞÞ2
! dB
Here, re, s, and ˆs represent the degraded, original and estimated exposure domain images, respectively. Figure 2a shows an image of a house. Figure 2b shows the image degraded by filmgrain noise. Figures 2c–f show images estimated by MWF, PF, ARUKF, and ISUKF.
58 Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
The ARUKF preserves sharp details (such as grass and tree leaves) akin to MWF, but is comparatively less noisy. The performance of PF is comparable to that of ARUKF. However, the image estimated by ISUKF is the best. The noise has been effectively filtered while simultaneously preserving fine details such as windowpanes, grass, and leaves. Note that qualitatively, the output of the proposed ISUKF method is strikingly good. This is also reflected in its high ISNR value. Next, we tested their performance on a large data set. For this purpose, we took a data set of 25 images that included faces, natural scenes, movie frames, textures, and satellite images, and compared the performance of ARUKF and ISUKF with MWF and PF. We performed the experiments with both moderate- and highnoise levels in several trials. Figure 3 gives the results in terms of mean ISNR values. From the plots, we observe that the proposed ISUKF outperforms other filters. ARUKF performance is better than MWFs and comparable to PF. The PF and ARUKF require the AR parameters of the
original image, which is a difficult proposition in real situations.
Figure 3. Performance comparison of filters on different images
Real degradation We now demonstrate performance on images degraded by real film-grain noise. Figure 4a shows a cropped portion of an image frame from the old classic movie The Testament of Dr. Mabuse. The film-grain noise shows up clearly on the white coat, the face, and so on. Figures 4b–e show the output of MWF, PF, and the ARUKF and ISUKF methods. We observed that ARUKF performs on par with PF, but the ISUKF method is the most effective even in real situations. The overall noise level is quite low. The folds on the coat are best brought out in the ISUKF’s output. Next, Figure 5a (next page) shows a face image with film-grain noise. Residual film-grain noise is visible in the MWF’s output (see Figure 5b). The PF result (see Figure 5c) is comparatively better. The ARUKF (see Figure 5d) is marginally better than PF. The ISUKF’s output (see Figure 5e) has all the sharp edges intact, and
over 20 Monte Carlo runs at: (a) moderate noise (s v2 5 0.05) and (b) high noise (s v2 5 0.15).
Figure 4. (a) Cropped portion of a frame from the movie The Testament of Dr. Mabuse. Output of (b) modified Wiener filter, (c) particle filter, (d) autoregressive unscented Kalman filter, and (e) importance sampling unscented Kalman filter. (Courtesy Criterion Collection.)
59 Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
Figure 5. (a) A face image with film-grain noise. Image estimated using (b) modified Wiener filter, (c) particle filter, (d) autoregressive unscented Kalman filter, and (e) importance sampling unscented Kalman filter.
Figure 6. (a) A building image. Output image obtained using a (b) modified Wiener filter, (c) particle filter, (d) autoregressive unscented Kalman filter, and (e) importance sampling unscented Kalman filter.
is least affected by grains. Finer details such as lips and eyebrows are well preserved. The image obtained using ISUKF is visually striking in quality and appears the most natural among them all. Finally, in Figure 6a we show another image with real film-grain noise. Figures 6b–e depict the output from MWF, PF, ARUKF, and ISUKF. The result obtained using ARUKF is comparable to PF. The ISUKF output is again the best. The edges are well preserved and the grainy appearance in the original degraded image (see Figure 6a) is greatly mitigated. Unlike the particle filter, in which we had to propagate about 200 samples throughout the filter in an attempt to faithfully represent the posterior distribution, for the ISUKF method we found that as few as 50 samples are sufficient to reliably estimate the first two moments of the non-Gaussian prior during prediction. At each pixel, the PF with S samples requires O(S) operations in weight calculations and O(S log S) comparisons in the resampling step. The ISUKF requires O(L) operations in the imporn3 n2 tance sampling step, 6a z 2a multiplications and additions for sigma point calculation, and O(2na + 1) operations for updating. Table 1 (on page 62) gives a summary of the filters’ computational requirements. The MWF is
simple and requires only O(MN log MN) operations. On a Pentium 4 PC with 256 Mbytes of RAM and for a 200 3 200 image, PF takes 148 seconds (for S 5 200). In comparison, the ISUKF (for L 5 100) and the ARUKF execute in 18 seconds and 14 seconds, respectively.
Conclusions In this article, we have discussed the application of the unscented Kalman filter for removal of film-grain noise in images. In addition to film-grain noise, old films and photographs also suffer from damage due to physical and chemical effects, scratches, stains, and digital drop-outs. Image inpainting, a companion topic to film-grain noise removal, refers to the process of reconstructing such damaged portions in images. One important extension of the UKF discussed here is for joint recovery of images from scratches and film-grain noise. The UKF can be combined with image inpainting techniques to achieve this goal. Given the huge legacy of old movies captured on photographic film, there is enormous scope for application of the techniques we have discussed in this article. MM
60 Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
Related Work Here we provide a brief summary of various methods reported to suppress film-grain noise. Andrews et al. expand the nonlinear observation model into a Taylor series and derived an approximate filter for recovering the original image.1 Using the film’s D 2 log E curve, and assuming Gaussian models for image and noise statistics, Hunt derives a maximum a posteriori (MAP) estimate of the restored image is derived.2 The resulting nonlinear equations are solved by steepest descent using fast transform computations. Naderi and Sawchuk proposed a linear filter that incorporates noise’s signal-dependent nature by adaptively changing the noise variance.3 Its performance is better than the conventional Wiener filter. Froehlich, Walkup, and Krile consider image estimation in signal-dependent, film-grain noise based on different priors.4 These authors employ numerical integration and solve higher-order polynomials to find the MAP estimate. Tavildar, Guptha, and Guptha model the image as Gaussian and use a signal-independent transformation to reduce the MAP solution’s complexity.5 Tekalp and Pavlovic propose a modified Wiener filter in the exposure domain, which is tailored specifically to tackle filmgrain noise.6 By incorporating sensor nonlinearity into the restoration procedure, they demonstrate significant improvements over the traditional wiener filter. Yan and Hatzinakos estimate the noise model parameters using higher-order statistics and propose a Wiener-type filter.7 They obtain filter coefficients by minimizing a cumulant criterion based on the extended Wiener-Hopf equations. Stefano, White, and Collis propose an undecimated wavelet-domain technique for reducing film-grain noise.8 Using a set of training images, the wavelet coefficients are thresholded for different noise levels by optimizing a cost function related to visual quality. Argenti, Torricelli, and Alparone perform signal-dependent film-grain noise reduction in both spatial and wavelet domains.9 The detail wavelet coefficients are adaptively shrunk using the local statistics derived from the noisy image. The wavelet filter is shown to outperform its spatial counterpart. Recently, Sadhar and Rajagopalan have demonstrated good results using a recursive spatial domain particle filter (PF)-based approach to tackle film-grain nonlinearity.10 The PF arrives at the conditional mean of the posterior by propagating samples and their associated weights within a Bayesian framework. However, the approach is computationally intensive and its performance is constrained by the assumption of a homogeneous autoregressive state model. Recursive estimation techniques have memory and implementation advantages and permit spatial adaptivity to be easily incorporated into the filter model.11 The dynamic state-space formulation of the Kalman filter can
readily incorporate an image’s space-varying local statistics into the estimation procedure. However, the Kalman filter is limited to linear state and measurement models. To deal with nonlinear systems, the extended Kalman filter (EKF) works by linearizing the state and measurement models. There have been recent successful works in 1D nonlinear filtering by employing the unscented Kalman filter (UKF) over the EKF.12 The UKF, proposed by Julier and Uhlmann, not only outperforms EKF in implementation ease and accuracy but is also more stable.12
References 1. H.C. Andrews and B.R. Hunt. Digital Image Restoration, PrenticeHall, 1977. 2. B.R. Hunt, ‘‘Bayesian Methods in Nonlinear Digital Image Restoration,’’ IEEE Trans. Computers, vol. 26, no. 3, 1977, pp. 219-229. 3. F. Naderi and A.A. Sawchuk, ‘‘Estimation of Images Degraded by FilmGrain Noise,’’ J. Applied Optics, vol. 20, no. 8, 1978, pp. 1228-1237. 4. G.K. Froehlich, J.F. Walkup, and T.F. Krile, ‘‘Estimation in SignalDependent Film-Grain Noise,’’ Applied Optics, vol. 20, no. 20, 1981, pp. 3619-3626. 5. A.S. Tavildar, H.M. Guptha, and S.N. Guptha, ‘‘Maximum A Posteriori Estimation in Presence of Film-Grain Noise,’’ Signal Processing, North-Holland, vol. 8, no. 3, 1985, pp. 363-368. 6. A.M. Tekalp and G. Pavlovic, ‘‘Restoration in the Presence of Multiplicative Noise with Application to Scanned Photographic Images,’’ IEEE Trans. Signal Processing, vol. 39, no. 9, 1991, pp. 2132-2136. 7. J.C.K. Yan and D. Hatzinakos, ‘‘Signal-Dependent Film Grain Noise Removal and Generation Based on Higher-Order Statistics,’’ Proc. IEEE Signal Processing Workshop on Higher-Order Statistics, 1997, pp. 77-82. 8. A.D. Stefano, P.R. White, and W.B. Collis, ‘‘Film Grain Reduction on Colour Images Using Undecimated Wavelet Transform,’’ Image and Vision Computing, vol. 22, no. 11, 2004, pp. 873-882. 9. F. Argenti, G. Torricelli, and L. Alparone, ‘‘MMSE Filtering of Generalised Signal-Dependent Noise in Spatial and Shift Invariant Wavelet Domains,’’ Signal Processing, vol. 86, no. 8, 2006, pp. 2056-2066. 10. S.I. Sadhar and A.N. Rajagopalan, ‘‘Image Recovery Under Nonlinear and Non-Gaussian Degradations,’’ J. Optical Society of America-A, vol. 22, no. 4, 2005, pp. 604-615. 11. E.S. Meinel, ‘‘Origins of Linear and Nonlinear Recursive Restoration Algorithms,’’ J. Optical Society of America–A, vol. 3, no. 6, 1986, pp. 787-799. 12. S.J. Julier and J.K. Uhlmann, ‘‘A New Extension of the Kalman Filter to Nonlinear Systems,’’ SPIE AeroSense Symp., 1997, pp. 182-193, http://www.cs.unc.edu/˜welch/kalman/media/ pdf/Julier1997_SPIE_KF.pdf.
61 Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
Table 1. Comparison of computational complexity at each pixel. Random Filter Particle filter Auto-regressive UKF Importance sampling UKF
Log and exponential evaluations 2S 2na + 1
Additions and multiplications
Inequality number comparisons generations O(S log S) 2S 0 0
O(S) n3a n2 z a z Oð2na z 1Þ 6 2
2L + 2na + 1 OðLÞ z
References 1. D. Zou et al., ‘‘Data Hiding in Film Grain,’’ Proc. 5th Int’l Workshop Digital Watermarking, LNCS 4283, Springer, 2006, pp. 197-211. 2. T.M. Moldovan, S. Roth, and M.J. Black, ‘‘Denois-
n3a n2 z a z Oð2na z 1Þ 6 2
0
L
11. F.C. Jeng and J.W. Woods, ‘‘Compound GaussMarkov Random Fields for Image Estimation,’’ IEEE Trans. Signal Processing, vol. 39, no. 3, 1991, pp. 683-697. 12. M. Ceccarelli, ‘‘Fast Edge-Preserving Picture Re-
ing Archival Films Using a Learned Bayesian
covery by Finite Markov Random Fields,’’ F. Roil
Model,’’ Proc. IEEE Int’l Conf. Image Processing,
and S. Vitulano, eds., Proc. Int’l Conf. Image Analysis
2006, pp. 2641-2644; http://www.cs.brown.edu/
and Processing, Springer Verlag, 2005, pp. 277-
people/roth/pubs/icip2006.pdf.
286.
3. H.C. Andrews and B.R. Hunt, Digital Image Restoration, Prentice-Hall, 1977.
13. S.Z. Li, Markov Random Field Modeling in Computer Vision, Springer Verlag, 1995.
4. A.M. Tekalp and G. Pavlovic, ‘‘Restoration in the
14. S. Geman and D. Geman, ‘‘Stochastic Relaxation,
Presence of Multiplicative Noise with Application
Gibbs Distributions, and the Bayesian Restoration
to Scanned Photographic Images,’’ IEEE Trans.
of Images,’’ IEEE Trans. Pattern Analysis and
Signal Processing, vol. 39, 1991, pp. 2132-2136.
Machine Intelligence., vol. 6, no. 11, 1984,
5. S.J. Julier and J.K. Uhlmann, ‘‘A New Extension of the Kalman Filter to Nonlinear Systems,’’ SPIE
pp. 721-741. 15. R. G. Aykroyd, ‘‘Bayesian Estimation for
AeroSense Symp., 1997, pp. 182-193; http://www. cs.unc.edu/˜welch/kalman/media/pdf/
Homogeneous and Inhomogeneous Gaussian
Julier1997_SPIE_KF.pdf.
Machine Intelligence, vol. 20, 1998, pp. 533-
6. R. van der Merwe et al., The Unscented Particle Filter, tech. Report CUED/F-INFENG/TR 380, Engineering Dept., Cambridge Univ., 2000. 7. S.J. Julier and J.K. Uhlmann, A General Method for Approximating Nonlinear Transformations of Probability Distributions, tech. report, RRG, Dept. of Engineering Science, Univ. of Oxford, 1996. 8. H. Kaufman et al., ‘‘Estimation and Identification of Two-Dimensional Images,’’ IEEE Trans. Automatic
Random Fields,’’ IEEE Trans. Pattern Analysis and 539. 16. D.J.C. Mackay, ‘‘Introduction to Monte Carlo Methods,’’ Learning in Graphical Models, M.I. Jordan, ed., NATO Science Series, Kluwer Academic Press, 1998, pp. 175–204. 17. S.I. Sadhar and A.N. Rajagopalan, ‘‘Image Recovery Under Nonlinear and Non-Gaussian Degradations,’’ J. Optical Society of America-A, vol. 22, 2005, pp. 604-615.
Control, vol. 28, no. 7, 1983, pp. 745-756.
IEEE MultiMedia
9. F.C. Jeng and J.W. Woods, ‘‘Inhomogeneous Gaussian Image Models for Estimation and
G.R.K. Sai Subrahmanyam
Restoration,’’ IEEE Trans. Acoustics Speech and
is a PhD candidate in the De-
Signal Proc., vol. 36, no. 8, 1988, pp. 1305-
partment of Electrical Engineer-
1312.
ing at the Indian Institute of
10. S.R. Kadaba, S.B. Gelfand, and R.L. Kashyap, ‘‘Recursive Estimation of Images Using NonGaussian Autoregressive Models,’’ IEEE Trans. Im-
Technology
Madras.
His
re-
search interests include image restoration, inpainting, and com-
age Processing, vol. 7, no. 10, 1998, pp. 1439-
puter vision. Subrahmanyam has an MS in electrical
1452.
engineering from the Indian Institute of Science,
62 Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
Bangalore.
Contact
him
at
saisubrahmanyam.
[email protected].
of the IEEE Transactions on Pattern Analysis and Machine Intelligence. Contact him at
[email protected]. ac.in.
A.N. Rajagopalan is an associate professor in the Depart-
Rangarajan Aravind is a fac-
ment of Electrical Engineering
ulty member in the Department
at
of
of Electrical Engineering at the
Technology Madras. His re-
Indian Institute of Technology
search interests include shape
Madras. Aravind has a PhD in
from
electrical engineering from the
the
X,
Indian
image
Institute
restoration,
super-resolution, and biometric recognition. Ra-
University of California, Santa
jagopalan has a PhD in electrical engineering
Barbara. His research interests include image and
from the Indian Institute of Technology Bombay.
video processing and compression. Contact him at
He is a Humboldt Fellow and an associate editor
[email protected].
April–June 2008 63 Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
