Importance Sampling-Based Unscented Kalman Filter for Film-Grain Noise Removal

12 Pages·2009·1.2 MB·English

Feature Article
In the density domain, an image affected by
Importance
film-grainnoisecanbemodeledasalogarith-
mic nonlinearity of the original scene cor-
rupted by additive white Gaussian noise.3,4
That is, r (m,n) 5 alog (s(m,n)) + b + v(m,n)
Sampling-Based d 10
where s(m,n) is the original scene, r (m,n) is
d
the degraded observation in the density
domain, and v(m, n) denotes additive white
Gaussian noise with zero mean and variance
Unscented
s2.Parametersaandbaretheslopeandoffset
v
derived from the film’s D 2 logE curve.
Alternatively, in the exposure domain, we
can write
Kalman Filter for
(cid:1) (cid:2)
r ðm,nÞ~sðm,nÞ 10vðm,nÞ=a ð1Þ
e
Film-Grain Noise
but the noise becomes multiplicative and
non-Gaussian. In Equation 1, r(m,n) is the
e
degraded exposure domain image and is
related to the density domain observation
Removal
image as r (m,n) 5 alog (r(m,n)) + b.
d 10 e
TheunscentedKalmanfilter(UKF)canincor-
porate multiplicative, non-Gaussian noise into
the estimation procedure. Demonstrating this
capabilityinthe2Ddomain,weproposeanauto-
G.R.K. Sai Subrahmanyam, A.N. Rajagopalan, and
regressive UKF for film-grain denoising. When
Rangarajan Aravind
modeled as an AR process, the image prior has
Indian Institute of Technology Madras
limitedcapabilitytohandleedges.Weproposea
novelmethodologythatreducesfilm-grainnoise
P
hotographic film is a widely used butpreservestheedgeinformationbyjudiciously
Photographicfilm
image-recording medium. Nega- incorporating a non-Gaussian prior within the
containsfilm-grain
tives contain tiny crystals of silver (UKF)frameworkthroughimportancesampling.
noisethattranslates
tomultiplicative, halide salts, which are light sensi- Weachievethisbyformulatingadiscontinuity-
non-Gaussian noise tive. When such a film is developed, these adaptive Markov random field (DAMRF) prior
inthe exposure crystals turn into tiny filaments of metallic suitedforrecursivepredictionandbyemploying
domain. Amethod silver conventionally called grains. These ran- importance sampling with a space-varying and
based onthe domlypatternedgrains,whichconstitutefilm- heavy-tailed proposal density to estimate the
unscentedKalman grain noise, are often visibly objectionable in DAMRFstatistics.
filter cansuppress photographic prints. While film grain can
thisnoisewhile enhance the feel of a motion picture,1 it State-space formulation and
simultaneously rendersthetaskofcompression(transmission) unscented transformation
preserving edge difficult because of the noise. For the multi-
TheUKFframeworkissimilartothatofthe
information. media industry, reducing film-grain noise is
Kalman filter. Like the Kalman filter, the UKF
important for scanning, digitization, and
isbasedonstate-spaceformulation butdiffers
efficient storage.2
in the way it represents and propagates
There is a well-known nonlinear relation-
Gaussianity through system dynamics. If the
ship between the incoming light intensity
image is modeled as a 2D autoregressive (AR)
(exposure) and the silver density deposited
process, the corresponding AR equation writ-
on the film. The domain in which the actual
teninstatetransitionform is
filmformationtakesplaceisreferredtoasthe
ddeenvseiltoypdeodmiaminagaendisthaveadiloambleainfoirnvwishuaiclhcothne- xðm,nÞ ~f(cid:3)xðm,n{1Þ,uðm,nÞ(cid:4)
sumptionisreferredtoastheexposuredomain. ~Fxðm,n{1Þzuðm,nÞ
52 1070-986X/08/$25.00G2008IEEE PublishedbytheIEEEComputerSociety
Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
Here,x(m,n)isthestatevectorwithdimension X ~xx(cid:1); wðmÞ~ l
nx. If we consider a first-order, causal, non- 0 0 ðnx zlÞ
symmetrichalfplane(NSHP)support,thennx X ~xx(cid:1)z(cid:1)pffiðffinffiffiffiffiffiffizffiffiffiffiffiffiffilffiffiÞffiffiPffiffiffi(cid:2),i~1,...,n ;
53,x(m,n)5[s(m,n),s(m,n21),s(m21,n)]T i x i x
andx(m,n21)5[s(m,n21),s(m21,n),s(m21, wðcÞ ~wðmÞz(cid:3)1{a2 zb (cid:4)
n21)]Twhere(m,n)denotesthecurrentpixel 0 0 UT UT
position,andsuperscriptTreferstotranspose. Xi~xx(cid:1){(cid:1)pðffiffinffiffiffixffiffiffizffiffiffiffiffiffiffilffiffiÞffiffiPffiffiffi(cid:2),i~nxz1,...,2nx;
i
The state transition matrix F contains the AR
1
coefficients of the image. Vector u 5 wðmÞ ~wðcÞ~ ,i~1,...,2n
(m,n) i i 2ðn zlÞ x ð2Þ
(u 0 0)T where u is state noise with x
(m,n) (m,n)
variance s2 and dimensionn (5 1).
u u
For the film-grain noise problem, the Here, l~a2UTðnxzkÞ{nx. Parameter aUT
measurement equation in exposure domain controls the spread of the sigma points
(seeEquation 1)isof theform around x¯ and is usually set to a value
between 0 and 1. The term b is used to
UT
incorporate prior knowledge about x. Note
yðm,nÞ ~h(cid:3)xðm,nÞ,vðm,nÞ(cid:4)~Hxðm,nÞ(cid:1)10vðm,nÞ=a(cid:2) that (cid:3)pffiðffinffiffiffixffiffiffizffiffiffiffiffiffiffilffiffiÞffiffiPffiffiffi(cid:4)i represents the ith row (or
column) of the n 3 n -matrix square root of
x x
(n + l)P.7 The addition and subtraction of
x
these vectors from the mean value of the
prior results in 2n symmetric sigma points.
x
wherey(m,n)5re(m,n)isascalarmeasurement Along with the prior mean X0 we have 2nx +
(with dimension ny 5 1), v(m, n) 5 v(m, n) is 1 sigma points. Subscript i on X refers to the
measurementnoiseofdimensionnv(51),and ith sigma point. The superscripts m and c on
H5[100].Forthefilm-grainnoiseproblem, the weights refer to their use in mean and
the measurement model h is nonlinear. covariance calculations, respectively. The
Hence, we resort to a nonlinear counterpart weight wðmÞ associated with the ith point
of theKalman filter,namely theUKF. i
satisfies P2nxwðmÞ~1. Simple matrix calcu-
The UKF is a straightforward application i~0 i
of the unscented transformation (UT) for lations can confirm that the first two mo-
recursive minimum mean-square-error esti- ments of these sigma points are exactly the
mation.5,6 The UT is founded on the notion same as that of the prior.
that with a fixed number of parameters, it’s We instantiate each sigma point through
easier to approximate a Gaussian distribu- the true nonlinearity to yield Yi 5 g(Xi),i 5
tion than to approximate an arbitrary non- 0,1,…,2nx. We estimate the mean and covari-
linear function or transformation.5,7 The UT ance of yas
is a deterministic sampling approach for
calculating the statistics of a random vari-
able that undergoes a nonlinear transforma- (cid:1)yy~ X2i~nx0wiðmÞYi
tion.
and
Considerpropagatingtherandomvariable
xthroughanonlinearfunctiong:Rnx?Rny to
generate y 5 g(x). In our problem, the Py ~ X2i~nx0wiðcÞðYi{(cid:1)yyÞðYi{(cid:1)yyÞT
function g is equivalent to f for the state
transformation, and h for the measurement
transformation. Assume x has mean x¯ and Positive semidefiniteness of the predicted
covariance P. To calculate the first two covariance matrix P is guaranteed by choos-
y
moments of y using the scaled UT, we ing k $ 0.6,7 These estimates of the moments
A
proceedasfollows:First,wedeterministically of y are accurate to the second-order (third- p
r
choose a set of 2nx + 1 weighted samples or orderforsymmetricpriors)oftheTaylorseries il–
J
sigmapointssothattheyexactlycapturethe expansion for any nonlinear function g(x), u
n
true mean and covariance of the prior because the prior sigma points completely e
2
random variable x. A selection scheme that capturethepriordistributionuptothesecond 0
satisfies this requirement,5,6 is as follows: moment.7 08
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 consists of sigma points
ðm,n{1Þ
Inthissection,weproposearecursivefilter, formed from the first n rows of Xa ,
x ðm,n{1Þ
basedontheUKF,thataccountsfortheimage whileXu isformedfromthen rowsthat
ðm,n{1Þ u
sensor nonlinearity intheobservation model.
immediately follow the first n rows in
The UKF, based on the UT, provides a x
Xa .
mathematically tractable way to propagate ðm,n{1Þ
Forstep4,weusethepredictedstatesigma
the first two moments even in the presence
points to determine the predicted mean and
of multiplicative and non-Gaussian noise. We
covariances:
achieve this by augmenting the state-random
vraanridaobmlewviathriatbhlees(staantedanapdpolybisnegrvaUtiToni)nnotihsee xx(cid:1)ðm,nÞ=ðm,n{1Þ ~ X2i~na0wiðmÞXxi,ðm,nÞ=ðm,n{1Þ
predictionofstateandobservation.5Weapply Pðm,nÞ=ðm,n{1Þ ~ X2i~na0wiðcÞ
the scaled UT sigma point selection
h i
scheme to the new augmented vector | Xxi,ðm,nÞ=ðm,n{1Þ {xx(cid:1)ðm,nÞ=ðm,n{1Þ
h iT
sxiaðomn,nÞn~ txoTðmc,naÞlcuuðmla,tneÞtvhðme,ncoÞrrespwointhdindgimaeung-- |hXxi,ðm,nÞ=ðm,n{1Þ {xx(cid:1)ðm,nÞ=ðm,n{1ÞiT
a
mented sigma point matrix Xa of dimen-
ðm,nÞ
sionn 3(2n +1).Heren 5n +n +n .(For
a a a x u v
example, for a first-order NSHP support where i represents the ith sigma point that is
nx~3,Xaðm,nÞ will be of size 5 3 11 and each theith columnof Xx .
sigma point isof dimension 5). ðm,nÞ=ðm,n{1Þ
Instep5,thesigmapointscorrespondingto
We now propose the auto-regressive UKF
themeasurementarepredictedbypropagating
(ARUKF), which sequentially estimates the in-
the (predicted) state sigma points and mea-
tensityoftheoriginalimageateachpixel(m,n)
surementnoisesigmapointsthroughthefilm-
inarasterscanorderfromlefttoright.Thefilter
grainnonlinearity intheexposure domain:
updatesthemeanandcovarianceoftheGauss-
ianapproximationtotheposteriordistribution (cid:1) (cid:2)
ofthestateaswedescribeinthefollowingsteps. Yðm,nÞ=ðm,n{1Þ~h Xxðm,nÞ=ðm,n{1Þ,Xvðm,n{1Þ
E[xW]e;Psta=rtEb[y(xin2itiax¯liz)(inxg2stax¯te)Tst]a.tistics x¯0 5 ~HXxðm,nÞ=ðm,n{1Þ:(cid:1)10:(cid:3)Xvðm,n{1Þ(cid:6)a(cid:4)
0 0 0 0 0 0
Step1augmentsstatemeanandcovariance
Here, Xv is formed from the last n rows
with noisestatistics: ðm,n{1Þ v
of Xa . Operations .* and .() represent
ðm,n{1Þ
h iT element-by-element operations (multiplica-
xx(cid:1)a ~ xx(cid:1)T 0 0
ðm,n{1Þ ðm,n{1Þ tion and raised power, respectively). This
2Pðm,n{1Þ 0 0 3 results in measurement sigma point matrix
Paðm,n{1Þ~664 0 s2u 0 775 Y(mIn,ns)/t(emp,n62,1w)eofesdtiimmeantesitohnenstya3tis(t2icnsar+eq1u).ired
0 0 s2
v forupdatingusingthepredictedsigmapoints:
For step 2, we calculate the augmented
sigma point matrix(seeEquation 2): (cid:1)yyðm,nÞ=ðm,n{1Þ ~ X2i~na0wiðmÞYi,ðm,nÞ=ðm,n{1Þ
Xaðm,n{1Þ ~ Pyy~X2i~na0wiðcÞhYi,ðm,nÞ=ðm,n{1Þ{(cid:1)yyðm,nÞ=ðm,n{1Þi
h qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffii
xx(cid:1)aðm,n{1Þxx(cid:1)aðm,n{1Þ+ ðnazlÞPaðm,n{1Þ |hYi,ðm,nÞ=ðm,n{1Þ {(cid:1)yyðm,nÞ=ðm,n{1ÞiT
a In step 3, the state sigma points and state Pxy~X2i~na0wiðcÞhXxi,ðm,nÞ=ðm,n{1Þ{xx(cid:1)ðm,nÞ=ðm,n{1Þi
Medi nAoRissetastiegmmaodpeolinfotsraprreepdricotpioagna:tedthroughthe |hYi,ðm,nÞ=ðm,n{1Þ {(cid:1)yyðm,nÞ=ðm,n{1ÞiT
ti
Mul Xxðm,nÞ=ðm,n{1Þ~f(cid:1)Xxðm,n{1Þ,Xuðm,n{1Þ(cid:2) where Yi, (m, n)/(m, n21) represents the ith
E sigma point, which is the ith column of
IEE ~FXxðm,n{1Þ zhXuðmT,n{1Þ 0T 0TiT 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 of the signal derivative exceeds a certain
filter: threshold.14 However, the resulting poten-
tial function becomes discontinuous. Li pro-
TheKalmangainKðm,nÞ~PxyP{yy1 posesacontinuousadaptiveregularizermodel
xx(cid:1)ðm,nÞ~xx(cid:1)ðm,nÞ=ðm,n{1Þ in which, unlike the line-field model, the
interaction decreases as the derivative magni-
(cid:1) (cid:2)
zKðm,nÞ yðm,nÞ{(cid:1)yyðm,nÞ=ðm,n{1Þ tude becomes larger and is completely turned
off at infinite magnitude of the signal deriva-
Pðm,nÞ ~Pðm,nÞ=ðm,n{1Þ{Kðm,nÞPyyKTðm,nÞ tive.13
FollowingLi,wemodelsasanon-Gaussian
Estimated pixel intensity at (m, n) isˆs(m,n)
MRFwith conditional PDFgiven by
5x¯ (1),thefirst component of x¯ .
(m,n) (m,n)
Pðsðm,nÞ=^ssðm{i,n{jÞÞ~
Non-Gaussian MRF prior
1 (cid:7) (cid:7) g2ðsðm,nÞ,^ssÞ(cid:8)(cid:8)
A homogeneous and linear AR model exp {clog 1z
Z c ð3Þ
imposesastrongconstraintonstateevolution,
which can hinder true transitions at image
where Z isa normalization constantand
edges.Otherresearchhasaddressedthisissue.
Forexample,theresearchofWoodsetal.uses
g2ðsðm,nÞ,^ssÞ~
spatiallyvarying2DARparameters, estimated
by windowing the observed image.8 Jeng and 1 X ðsðm,nÞ{^ssðm{i,n{jÞÞ2
Woods propose inhomogeneous AR image r2 ði2 zj2Þ
models using local statistics.9 Kadaba et al. ðm,nÞði,jÞ[R
observe that Bernouli-Laplacian and Cauchy Here,(i,j)MR5{(i,j)I(1#i#M ,2M #j#
1 1
distribution are a better fit than the Gaussian
M ) < (0, 1# i #M )},and M is theorder of
for error residuals.10 Jeng and Woods use a 1 1 1
the NSHP support. Figure 1a (next page)
compound Gauss MRF to model images.11
illustrates NSHP support. Parameter r2
Ceccarelliproposededge-preservingMRFs.12 ðm,nÞ
controlstheallowablevariationofthecurrent
We base our approach to incorporating a
pixel with respect to its neighboring values.15
non-Gaussian prior within our recursive esti-
Here, thepotentialfunction that accounts for
mation scheme on the fact that statistical
the irregularities in the solution is given by
dependenceincorporatedusinganMRFmodel
(cid:1) (cid:2)
provides more flexibility in incorporating clog 1z g2 .TheplotinFigure1bshowsthe
c
contextual constraints through a suitably
degree of penalty dictated by the potential
derived prior. An MRF possesses a Markovian
function. For this MRF model, the smoothing
property—thatis,thevalueofapixeldepends
strengthoftheregularizerincreasesmonoton-
onlyonthevaluesofitsneighboringpixels.13
ically as g increases within the convex band
It’s important to realize that at discontinu- B ~(cid:3){pfficffiffi,pfficffiffi(cid:4).Outsidetheband,smoothing
c
ities or edges, the notion of neighborhood
decreasesandbecomeszeroasgR‘.Because
dependencyshouldberelaxed.Incorporating
this feature preserves image discontinuities,
the smoothness constraint while simulta-
it’scalled aDAMRFmodel.13
neouslypreservingtheedgesisachallenging
task.
For simplicity and analytical tractability, it Importance sampling
isnormalpracticetomodeltheoriginalimage Toincorporatethenon-Gaussianpriorgiven
s by a homogeneous Gaussian MRF. The issue bytheDAMRFmodel(seeEquation3)intothe
with a homogeneous Gaussian MRF is that it UKFframework,weneedtoestimatethemean
imposes smoothness constraints everywhere, and covariance of this prior. Because it is
which inevitably leads to oversmoothing edg- analytically impossible to compute these mo-
A
es. A better way to handle the situation is to ments, we resort to a Monte Carlo method p
r
use a non-Gaussian conditional probability known as importance sampling.16 This tech- il–
J
density function (PDF) to model the original niquedeterminesatargetPDF’sestimatesusing u
n
image.GemanandGemanproposetheuseof thesamplesfromaneasy-to-sampleapproxima- e
2
linefieldsinwhichthesmoothnessconstraint tion.ConsideraPDFp(z)thatisknownuptoa 0
0
isswitchedoffatpointswherethemagnitude multiplicative constant but is difficult to esti- 8
55
Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
Figure1.(a) ð ð (cid:7)pðzÞ(cid:8)
Nonsymmetrichalf- Ep½gðzÞ(cid:2)~ gðzÞpðzÞdz~ gðzÞ qðzÞ qðzÞdz
planesupportata
pixel(m,n).(b)The (cid:9) (cid:7)pðzÞ(cid:8)(cid:10) 1XL (cid:1) (cid:2) p(cid:3)zðlÞ(cid:4)!
~E gðzÞ & g zðlÞ
degreeofpenaltywith q qðzÞ L qðzðlÞÞ
l~1
distanceimposedby
thediscontinuity- Here, z(l) represents samples drawn from q(z).
adaptiveMarkov Thebasicideabehindimportancesamplingis
randomfieldmodel. to choose a distribution that emphasizes
(c)Importance certain values of the input random variables
sampling:pisthe thathavemoreimpactontheparameterbeing
targetprobability estimatedthanothersbysamplingthemmore
densityfunctionwhose frequently.
momentsaretobe Our aim is to determine the estimates
estimated,whileqis under the PDF p, using samples from q.
the(Cauchy)
Because it’s difficult to draw samples from
samplerdensity. thetargetPDFp,wedrawLsamples,(cid:11)zðlÞ(cid:12)L
l~1
from the sampler PDF q—that is, instead of
choosing samples from a uniform distribu-
tion,theyarenowchosenfromadistribution
qthatconcentratesonthosepointswherethe
functionpislarge.Whenweusesamplesfrom
q to determine any estimates under p, in the
regions where q is greater than p, the esti-
mates are over-represented. In the regions
where q is less than p, they are under-
represented. To account for this, we use
correction weights
p(cid:3)zðlÞ(cid:4)
wl~
qðzðlÞÞ
in determining the estimates under p. For
example,tofindthemeanofthedistribution
pweuse
PL wlzðlÞ
^m ~ l~1
p PL wl
l~1
As L R ‘, the estimate mˆ tends to the actual
p
mate its moments. From the functional form, meanvalueofp.Inourproblem,theDAMRF
we can estimate its support. Consider a distri- distribution corresponds to p. We choose the
butionq(z),areasonableapproximationtop(z), Cauchydistributionfortheimportancefunc-
which is also known up to a multiplicative tion q, as it yields better estimates (due to its
constant, is easy to sample, and the(non-zero) heavy tailed nature) compared to the Gauss-
supportofq(z)includesthesupportofp(z).Such iansampler.16
a density q(z) is called the sampler density.
a Figure1c shows a typical plot that illustrates Edge-preserving UKF
di thetargetPDFpandsamplerdistributionq. We now extend the basic structure of the
e
M Mathematically, we can explain theprinci- UKF to include a non-Gaussian prior. In the
ulti ple of importance sampling as follows. Sup- proposed strategy, we use knowledge of the
M posethedensityq(z)roughlyapproximatesthe conditional PDF to capture the image’s local
E density of interest p(z), then the expected physical characteristics. Doing so implicitly
E
E valueofafunctiongofarandomvariablezis generalizes the state transition equation and
I
56
Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
doesn’trestrictittobeing linear. Webaseour The mean mˆp and variance s^2p of p are
proposed extension on the observation that
computed as
the Kalman filter variants differ only in state
vector propagation but have same update PL wlzðlÞ
equations.5 m^p ~ Pl~L1 wl
l~1
Specifically, we formulate the state condi-
tional density of the pixel to be estimated and
(given its neighbors) based on the DAMRF
model discussed previously. We estimate the PL wl(cid:1)zðlÞ {m^ (cid:2)2
mean and variance of the non-Gaussian s^2 ~ l~1 p
conditional prior, which constitutes the pre- p PLl~1wl
dicted statistics of the state, by importance
Instep3,wedirectlyusetheseestimatesto
sampling. On the basis of these and known
predict the one-step-ahead sigma points.
measurement noise statistics, we determine
These estimates for the mean and covariance
predicted sigma points and propagate them
are augmented to deal with the nonadditive
through the nonlinearity. The approach takes
noise,as thefollowing shows:
the sigma points to the update stage of the
UKFanddetermines thefinal imageestimate.
The steps involved in the proposed impor- xx(cid:1)ðm,nÞ=ðm,n{1Þ~ ^mp,
tancesamplingUKF(ISUKF)methodfollows. Pðm,nÞ=ðm,n{1Þ~s^2p,
Forstep1,ateachpixel,constructthestate
h iT
conditional PDF using the past pixels in the xx(cid:1)a ~ xx(cid:1)T 0
ðm,nÞ=ðm,n{1Þ ðm,nÞ=ðm,n{1Þ
NSHPsupport,andthevaluesofr2 andcin
ðm,nÞ h iT
theDAMRF model.That is, ~ ^m 0 ,and
p
Pðsðm,nÞ=^ssðm{i,m{jÞÞ "Pðm,nÞ=ðm,n{1Þ 0 #
Pa ~
~exp(cid:7){clog(cid:7)1z g2ðsðm,nÞ,^ssÞ(cid:8)(cid:8) ðm,nÞ=ðm,n{1Þ 0 s2v
c 2s^2 0 3
~ p
4 5
If M 5 1, then the neighbors considered in 0 s2
1 v
the first-order NSHP support are ˆs(m,n 2 1),
ˆs(m21,n),ˆs(m21,n21),andˆs(m21,n+1),
Because we base prediction directly on the
which already are estimated. Because the
priorPDF,weneednotaugmentwiththestate
parameter r2ðm,nÞ of the DAMRF model charac- noise.
terizes the local dependency, we set it to Weuse these augmented meanand covari-
P(m,n21), thecovariance of thepreviousstate. ancestodeterminethesigmapointmatrix(see
Forstep2,obtainthemeanandcovariance Equation 2):
of the previous PDF using importance sam-
pling. Draw samples (cid:11)zðlÞ(cid:12)Ll~1 from a Cauchy Xaðm,nÞ=ðm,n{1Þ~
sampler. Cauchy distributed samples with h
xx(cid:1)a
locationk andscales canbegeneratedusing ðm,nÞ=ðm,n{1Þ
1 1
k + s tan(p(u 2 0.5)), where u is uniformly qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffii
1 1 xx(cid:1)a + ðn zlÞPa
distributed from 0 to 1. To facilitate a close ðm,nÞ=ðm,n{1Þ a ðm,nÞ=ðm,n{1Þ
match with the DAMRF PDF p at each pixel,
the location of the sampler density q(z) is
varied as mq~ 14ð^ssðm,n{1Þz^ssðm{1,nÞz Then we propagate each of the 2na + 1
^ssðm{1,n{1Þz^ssðm{1,nz1ÞÞ and the sigma points corresponding to state and
scale is chosen so that it is large enough to measurement noise through the film-grain A
p
includethesupportofpinq.Thesamples are nonlinearity. r
weightedby correctionweights (cid:1) (cid:2) il–Ju
p(cid:3)zðlÞ(cid:4) Yðm,nÞ=m,n{1Þ ~h Xxðm,nÞ=m,n{1Þ,Xvðm,n{1Þ ne
wl~ qðzðlÞÞ ~Xxðm,nÞ=m,n{1Þ:(cid:1)10:(cid:3)Xvðm,n{1Þ(cid:6)a(cid:4) 200
8
57
Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
Figure2.(a)Imageofa where Xxðm,nÞ=m,n{1Þ and Xvðm,n{1Þ are formed chose a 5 5 for all experiments. We set the
house.(b)Degraded fromthefirstnxrowsofXaðm,nÞ=m,n{1Þandthenv parametersasaUT51,bUT50,andk51.The
image(s2v ~0:05).(c) rows following nx, respectively. Here .* and.() UKF is robust to variations in these parame-
Resultusinga are element-by-element operations. ters.
modifiedWienerfilter Finally,weestimatethestatisticsofy(from
(improvementin thesigmapoints)andupdatefollowingsteps6 Simulated degradation
signal-to-noiseratio5 For a quantitative comparison, we use the
and7asintheauto-regressiveUKF.Wederive
2.51dB),(d)particle improvement in signal-to-noise ratio (ISNR),
the estimated pixel intensity by ˆs(m,n) 5
filter(ISNR5 definedas
x¯ . Thus,ˆsyields theestimatedimage.
(m,n)
3.48dB),(e)auto-
regressiveunscented Experimental results P ðr ðm,nÞ{sðm,nÞÞ2!
Kalmanfilter(ISNR5 ISNR~10log m,n e dB
3.40dB),and(f) Here,wepresentresultsobtainedusingthe 10 Pm,nð^ssðm,nÞ{sðm,nÞÞ2
proposed ARUKF and ISUKF methods and
importancesampling
compare their performance with the stand-
unscentedKalman
ard modified Wiener filter (MWF)4 and with Here, r, s, and ˆs represent the degraded,
filter(ISNR5 e
therecentlyproposedparticlefilter(PF)meth- original and estimated exposure domain im-
4.55dB).
od.17 ages,respectively.
The c and r(m,n) parameters influence the Figure2a shows an image of a house.
DAMRFmodel.Wehavefoundexperimentally Figure2b shows the image degraded by film-
thatc51.8yieldsthebestperformance.From grainnoise.Figures 2c–fshowimagesestimat-
the D 2 logE characteristics of the film, we edby 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 ARUKFpreserves sharp details (suchas originalimage,whichisadifficultproposition Figure3.Performance
grass and tree leaves) akin to MWF, but is inrealsituations. comparisonoffilters
comparatively less noisy. The performance of ondifferentimages
PF is comparable to that of ARUKF. However, Real degradation over20MonteCarlo
theimage estimatedbyISUKFisthebest.The We now demonstrate performance on im- runsat:(a)moderate
noisehasbeeneffectivelyfilteredwhilesimul- ages degraded by real film-grain noise. Fig- noise(sv250.05)and
taneously preserving fine details such as ure 4a shows a cropped portion of an image (b)highnoise(sv2
windowpanes, grass, and leaves. Note that framefromtheoldclassicmovieTheTestament 50.15).
qualitatively, the output of the proposed of Dr. Mabuse. The film-grain noise shows up
ISUKF method is strikingly good. This is also clearlyonthewhitecoat,theface,andsoon.
reflectedinits highISNR value. Figures 4b–eshowtheoutputofMWF,PF,and
Figure4.(a)Cropped
Next,wetestedtheirperformanceonalarge the ARUKF and ISUKF methods. We observed
portionofaframe
dataset.Forthispurpose,wetookadatasetof that ARUKF performs on par with PF, but the
fromthemovieThe
25 images that included faces, natural scenes, ISUKF method is the most effective even in
TestamentofDr.
movie frames, textures, and satellite images, real situations. The overall noise level is quite Mabuse.Outputof(b)
andcomparedtheperformanceofARUKFand low.Thefoldsonthecoatarebestbroughtout modifiedWienerfilter,
ISUKF with MWF and PF. We performed the intheISUKF’s output. (c)particlefilter,(d)
experiments with both moderate- and high- Next, Figure5a (next page) shows a face autoregressive
noise levels in several trials. Figure3 gives the imagewithfilm-grainnoise.Residualfilm-grain unscentedKalman
resultsintermsofmeanISNRvalues.Fromthe noise is visible in the MWF’s output (see filter,and(e)
plots, we observe that the proposed ISUKF Figure5b). The PF result (see Figure5c) is importancesampling
outperforms other filters. ARUKF performance comparativelybetter.TheARUKF(seeFigure5d) unscentedKalman
isbetterthanMWFsandcomparabletoPF.The ismarginallybetterthanPF.TheISUKF’soutput filter.(Courtesy
PFandARUKFrequiretheARparametersofthe (seeFigure5e)hasallthesharpedgesintact,and CriterionCollection.)
59
Authorized licensed use limited to: EPFL LAUSANNE. Downloaded on May 18, 2009 at 04:10 from IEEE Xplore. Restrictions apply.
Figure5.(a)Aface isleastaffectedbygrains.Finerdetailssuchaslips simple and requires only O(MN log MN)
imagewithfilm-grain and eyebrows are well preserved. The image operations. On a Pentium 4 PC with 256
noise.Imageestimated obtainedusingISUKFisvisuallystrikinginquality MbytesofRAMandfora2003200image,PF
using(b)modified andappearsthemostnaturalamongthemall. takes148seconds(forS5200).Incomparison,
Wienerfilter,(c) Finally, in Figure6a we show another theISUKF(forL5100)andtheARUKFexecute
particlefilter,(d)auto- image with real film-grain noise. Figures6b–e in18secondsand14seconds,respectively.
regressiveunscented depicttheoutputfromMWF,PF,ARUKF,and
Kalmanfilter,and(e) ISUKF. The result obtained using ARUKF is
Conclusions
importancesampling comparable to PF. The ISUKF output is again
In this article, we have discussed the
unscented thebest. Theedges are well preserved and the
application of the unscented Kalman filter
Kalmanfilter. grainy appearance in the original degraded
for removal of film-grain noise in images. In
image (see Figure6a) is greatly mitigated.
addition to film-grain noise, old films and
Unlike the particle filter, in which we had to
photographs also suffer from damage due to
propagate about 200 samples throughout the
physical and chemical effects, scratches,
filter in an attempt to faithfully represent the
posteriordistribution,fortheISUKFmethodwe stains, and digital drop-outs. Image inpaint-
foundthatasfewas50samplesaresufficientto ing, a companion topic to film-grain noise
Figure6.(a)Abuilding reliably estimate thefirst twomomentsof the removal,referstotheprocessofreconstructing
image.Outputimage non-Gaussian priorduringprediction.Ateach suchdamaged portions inimages.
obtainedusinga(b) pixel, the PF with S samples requires O(S) One important extension of the UKF dis-
modifiedWienerfilter, operations in weight calculations and O(S log cussed here is for joint recovery of images
(c)particlefilter,(d) S) comparisons in the resampling step. The fromscratchesandfilm-grainnoise.TheUKF
autoregressive ISUKF requires O(L) operations in the impor- can be combined with image inpainting
unscentedKalman tance sampling step, n3a z n2a multiplications techniques to achieve this goal. Given the
6 2
filter,and(e) andadditionsforsigmapointcalculation,and huge legacy of old movies captured on
importancesampling O(2n +1)operationsforupdating.Table1(on photographic film, there is enormous scope
a
unscented page 62) gives a summary of the filters’ for application of the techniques we have
Kalmanfilter. computational requirements. The MWF is discussedinthisarticle. 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 readily incorporate an image’s space-varying local statistics
reportedtosuppressfilm-grainnoise.Andrewsetal.expand intotheestimationprocedure.However,theKalmanfilteris
the nonlinear observation model into a Taylor series and limited to linear state and measurement models. To deal
derived an approximate filter for recovering the original with nonlinear systems, the extended Kalman filter (EKF)
image.1 Using the film’s D 2 log E curve, and assuming works by linearizing the state and measurement models.
Gaussianmodelsforimageandnoisestatistics,Huntderives There have been recent successful works in 1D nonlinear
a maximum a posteriori (MAP) estimate of the restored filtering by employing the unscented Kalman filter (UKF)
image is derived.2 The resulting nonlinear equations are over the EKF.12 The UKF, proposed by Julier and Uhlmann,
solved by steepest descent using fast transform computa- not only outperforms EKF in implementation ease and
tions. Naderi and Sawchuk proposed a linear filter that accuracybut isalsomore stable.12
incorporates noise’s signal-dependent nature by adaptively
changingthenoisevariance.3Itsperformanceisbetterthan References
the conventional Wiener filter. Froehlich, Walkup, and Krile 1. H.C.AndrewsandB.R.Hunt.DigitalImageRestoration,Prentice-
consider image estimation in signal-dependent, film-grain Hall,1977.
noise based on different priors.4 These authors employ 2. B.R.Hunt,‘‘BayesianMethodsinNonlinearDigitalImage
numericalintegrationandsolvehigher-orderpolynomialsto Restoration,’’IEEETrans.Computers,vol.26,no.3,1977,
findtheMAPestimate.Tavildar,Guptha,andGupthamodel pp.219-229.
the image as Gaussian and use a signal-independent 3. F.NaderiandA.A.Sawchuk,‘‘EstimationofImagesDegradedbyFilm-
transformation to reduce the MAP solution’s complexity.5 GrainNoise,’’J.AppliedOptics,vol.20,no.8,1978,pp.1228-1237.
Tekalp and Pavlovic propose a modified Wiener filter in the 4. G.K.Froehlich,J.F.Walkup,andT.F.Krile,‘‘EstimationinSignal-
exposuredomain,whichistailoredspecificallytotacklefilm- DependentFilm-GrainNoise,’’AppliedOptics,vol.20,no.20,
grain noise.6 By incorporating sensor nonlinearity into the 1981,pp.3619-3626.
restoration procedure, they demonstrate significant im- 5. A.S.Tavildar,H.M.Guptha,andS.N.Guptha,‘‘MaximumA
provements over the traditional wiener filter. Yan and PosterioriEstimationinPresenceofFilm-GrainNoise,’’Signal
Hatzinakos estimate the noise model parameters using Processing,North-Holland,vol.8,no.3,1985,pp.363-368.
higher-order statistics and propose a Wiener-type filter.7 6. A.M.TekalpandG.Pavlovic,‘‘RestorationinthePresenceof
They obtain filter coefficients by minimizing a cumulant MultiplicativeNoisewithApplicationtoScannedPhotographic
criterion basedon theextendedWiener-Hopf equations. Images,’’IEEETrans.SignalProcessing,vol.39,no.9,1991,
Stefano, White, and Collis propose an undecimated pp.2132-2136.
wavelet-domain technique for reducing film-grain noise.8 7. J.C.K.YanandD.Hatzinakos,‘‘Signal-DependentFilmGrain
Using a set of training images, the wavelet coefficients are NoiseRemovalandGenerationBasedonHigher-OrderStatis-
thresholded for different noise levels by optimizing a cost tics,’’Proc.IEEESignalProcessingWorkshoponHigher-Order
function related to visual quality. Argenti, Torricelli, and Statistics,1997,pp.77-82.
Alparone perform signal-dependent film-grain noise reduc- 8. A.D.Stefano,P.R.White,andW.B.Collis,‘‘FilmGrainReduction
tioninbothspatialandwaveletdomains.9Thedetailwavelet onColourImagesUsingUndecimatedWaveletTransform,’’
coefficients are adaptively shrunk using the local statistics ImageandVisionComputing,vol.22,no.11,2004,pp.873-882.
derivedfromthenoisyimage.Thewaveletfilterisshownto 9. F.Argenti,G.Torricelli,andL.Alparone,‘‘MMSEFilteringof
outperform its spatial counterpart. Recently, Sadhar and GeneralisedSignal-DependentNoiseinSpatialandShift
Rajagopalan have demonstrated good results using a InvariantWaveletDomains,’’SignalProcessing,vol.86,no.8,
recursive spatial domain particle filter (PF)-based approach 2006,pp.2056-2066.
to tackle film-grain nonlinearity.10 The PF arrives at the 10. S.I.SadharandA.N.Rajagopalan,‘‘ImageRecoveryUnder
conditional mean of the posterior by propagating samples NonlinearandNon-GaussianDegradations,’’J.OpticalSocietyof
and their associated weights within a Bayesian framework. America-A,vol.22,no.4,2005,pp.604-615.
However, the approach is computationally intensive and its 11. E.S.Meinel,‘‘OriginsofLinearandNonlinearRecursive
performance is constrained by the assumption of a homo- RestorationAlgorithms,’’J.OpticalSocietyofAmerica–A,vol.3,
geneous autoregressive statemodel. no.6,1986,pp.787-799.
Recursive estimation techniques have memory and 12. S.J.JulierandJ.K.Uhlmann,‘‘ANewExtensionoftheKalman
implementation advantages and permit spatial adaptivity FiltertoNonlinearSystems,’’SPIEAeroSenseSymp.,1997,
to be easily incorporated into the filter model.11 The pp.182-193,http://www.cs.unc.edu/ welch/kalman/media/
˜
dynamic state-space formulation of the Kalman filter can 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.