AFastApproximationoftheBilateralFilter
usingaSignalProcessingApproach
SylvainParisandFr′edoDurand
MassachusettsInstituteofTechnology
ComputerScienceandArtificialIntelligenceLaboratory
Abstract.Thebilateralfilterisanonlinearfilterthatsmoothesasignal
whilepreservingstrongedges.Ithasdemonstratedgreate?ectivenessfor
avarietyofproblemsincomputervisionandcomputergraphics,anda
fastversionhasbeenproposed.Unfortunately,littleisknownaboutthe
accuracyofsuchacceleration.Inthispaper,weproposeanewsignal-
processinganalysisofthebilateralfilter,whichcomplementstherecent
studiesthatanalyzeditasaPDEorasarobuststatisticsestimator.
Importantly,thissignal-processingperspectiveallowsustodevelopa
novelbilateralfilteringaccelerationusingadownsamplinginspaceand
intensity.Thisa?ordsaprincipledexpressionoftheaccuracyintermsof
bandwidthandsampling.Thekeytoouranalysisistoexpressthefilter
inahigher-dimensionalspacewherethesignalintensityisaddedtothe
originaldomaindimensions.Thebilateralfiltercanthenbeexpressed
assimplelinearconvolutionsinthisaugmentedspacefollowedbytwo
simplenonlinearities.Thisallowsustoderivesimplecriteriafordown-
samplingthekeyoperationsandtoachieveimportantaccelerationofthe
bilateralfilter.Weshowthat,forthesamerunningtime,ourmethodis
significantlymoreaccuratethanpreviousaccelerationtechniques.
1Introduction
ThebilateralfilterisanonlinearfilterproposedbyTomasiandManduchito
smoothimages[1].Ithasbeenadoptedforseveralapplicationssuchastexture
removal[2],dynamicrangecompression[3],andphotographenhancement[4,5].
Ithasalsobeadaptedtootherdomainssuchasmeshfairing[6,7],volumetric
denoising[8]andexposurecorrectionofvideos[9].Thislargesuccessstemsfrom
severalorigins.First,itsformulationandimplementationaresimple:apixelis
simplyreplacedbyaweightedmeanofitsneighbors.Anditiseasytoadaptto
agivencontextaslongasadistancecanbecomputedbetweentwopixelvalues
(e.g.distancebetweenhairorientationsin[10]).Thebilateralfilterisalsonon-
iterative,i.e.itachievessatisfyingresultswithonlyasinglepass.Thismakes
thefilter’sparametersrelativelyintuitivesincetheiractiondoesnotdependon
thecumulatede?ectsofseveraliterations.
Ontheotherhand,thebilateralfilterisnonlinearanditsevaluationiscom-
putationallyexpensivesincetraditionalacceleration,suchasperformingconvolu-
tionafteranFFT,arenotapplicable.Elad[11]proposesanaccelerationmethod
usingGauss-Seideliterations,butitonlyapplieswhenmultipleiterationsofthe
2
filterarerequired.DurandandDorsey[3]describealinearizedversionofthe
filterthatachievesdramaticspeed-upsbydownsamplingthedata,achieving
runningtimesunderonesecond.Unfortunately,thistechniqueisnotgrounded
onfirmtheoreticalfoundations,anditisdi?culttoevaluatetheaccuracythat
issacrificed.Inthispaper,webuildonthisworkbutweinterpretthebilateral
filterintermsofsignalprocessinginahigher-dimensionalspace.Thisallowsus
toderiveanimprovedaccelerationschemethatyieldsequivalentrunningtimes
butdramaticallyimprovesnumericalaccuracy.
ContributionsThispaperintroducesthefollowingcontributions:
–Aninterpretationofthebilateralfilterinasignalprocessingframework.
Usingahigherdimensionalspace,weformulatethebilateralfilterasacon-
volutionfollowedbysimplenonlinearities.
–Usingthishigherdimensionalspace,wedemonstratethattheconvolution
computationcanbedownsampledwithoutsignificantlyimpactingtheresult
accuracy.Thisapproximationtechniqueenablesaspeed-upofseveralorders
ofmagnitudewhilecontrollingtheerrorinduced.
2RelatedWork
ThebilateralfilterwasfirstintroducedbySmithandBradyunderthename
“SUSAN”[12].ItwasrediscoveredlaterbyTomasiandManduchi[1]whocalled
itthe“bilateralfilter”whichisnowthemostcommonlyusedname.Thefilter
replaceseachpixelbyaweightedaverageofitsneighbors.Theweightassigned
toeachneighbordecreaseswithboththedistanceintheimageplane(thespatial
domainS)andthedistanceontheintensityaxis(therangedomainR).Using
aGaussianG
σ
asadecreasingfunction,andconsideringagrey-levelimageI,
theresultI
bf
ofthebilateralfilterisdefinedby:
I
bf
p
=
1
W
bf
p
summationdisplay
q∈S
G
σ
s
(||p?q||)G
σ
r
(|I
p
?I
q
|)I
q
(1a)
withW
bf
p
=
summationdisplay
q∈S
G
σ
s
(||p?q||)G
σ
r
(|I
p
?I
q
|)(1b)
Theparameterσ
s
definesthesizeofthespatialneighborhoodusedtofiltera
pixel,andσ
r
controlshowmuchanadjacentpixelisdownweightedbecauseof
theintensitydi?erence.W
bf
normalizesthesumoftheweights.
Barash[13]showsthatthetwoweightfunctionsareactuallyequivalentto
asingleweightfunctionbasedonadistancedefinedonS×R.Usingthisap-
proach,herelatesthebilateralfiltertoadaptivesmoothing.Ourworkfollowsa
similarideaandalsousesS×Rtodescribebilateralfiltering.Ourformulationis
nonethelesssignificantlydi?erentbecausewenotonlyusethehigher-dimensional
spaceforthedefinitionofadistance,butwealsouseconvolutioninthisspace.
Elad[11]demonstratesthatthebilateralfilterissimilartotheJacobialgorithm,
withthespecificitythatitaccountsforalargerneighborhoodinsteadoftheclos-
estadjacentpixelsusuallyconsidered.Buadesetal.[14]exposeanasymptotic
analysisoftheYaroslavskyfilter[15]whichisaspecialcaseofthebilateralfil-
terwithastepfunctionasspatialweight.Theyprovethatasymptotically,the
3
YaroslavskyfilterbehavesasthePerona-Malikfilter,i.e.italternatesbetween
smoothingandshockformationdependingonthegradientintensity.Durandand
Dorsey[3]casttheirstudyintotherobuststatisticsframework[16,17].They
showthatthebilateralfilterisaw-estimator[17](p.116).Thisexplainsthe
roleoftherangeweightintermsofsensitivitytooutliers.Theyalsopointout
thatthebilateralfiltercanbeseenasanextensionofthePerona-Malikfilter
usingalargerneighborhood.Mr′azeketal.[18]relatebilateralfilteringtoalarge
familyofnonlinearfilters.Fromasingleequation,theyexpressfilterssuchas
anisotropicdi?usionandstatisticalestimatorsbyvaryingtheneighborhoodsize
andtheinvolvedfunctions.Themaindi?erencebetweenourstudyandexisting
workisthatthepreviousapproacheslinkbilateralfilteringtoanothernonlinear
filterbasedonPDEsorstatisticswhereaswecastourstudyintoasignalprocess-
ingframework.Wedemonstratethatthebilateralfiltercanbemainlycomputed
withlinearoperations,thenonlinearitiesbeinggroupedinafinalstep.
Severalarticles[11,14,19]improvethebilateralfilter.Theysharethesame
idea:Byexploitingthelocal“slope”oftheimageintensity,itispossibletobetter
representthelocalshapeofthesignal.Thus,theydefineamodifiedfilterthat
betterpreservetheimagecharacteristicse.g.theyavoidtheformationofshocks.
Wehavenotexploredthisdirectionsincetheformulationbecomessignificantly
morecomplex.Itishoweveraninterestingavenueforfuturework.
Theworkmostrelatedtooursarethespeed-uptechniquesproposedby
Elad[11]andDurandandDorsey[3].Elad[11]usesGauss-Seideliterationsto
acceleratetheconvergenceofiterativefiltering.Unfortunately,noresultsare
shown–andthistechniqueisonlyusefulwhenthefilterisiteratedtoreachthe
stablepoint,whichisnotitsstandarduseofthebilateralfilter(oneiteration
oronlyafew).DurandandDorsey[3]linearizethebilateralfilterandpropose
adownsamplingschemetoacceleratethecomputationdowntofewsecondsor
less.However,notheoreticalstudyisproposed,andtheaccuracyoftheapprox-
imationisunclear.Incomparison,webaseourtechniqueonsignalprocessing
groundswhichhelpustodefineanewandmeaningfulnumericalscheme.Our
algorithmperformslow-passfilteringinahigher-dimensionalspacethanDurand
andDorsey’s[3].Thecostofahigher-dimensionalconvolutioniso?setbythe
accuracygain,whichyieldsbetterperformanceforthesameaccuracy.
3SignalProcessingApproach
Wedecomposethebilateralfilterintoaconvolutionfollowedbytwononlineari-
ties.Tocastthefilterasaconvolution,wedefineahomogeneousintensitythat
willallowustoobtainthenormalizationtermW
bf
p
asanhomogeneouscompo-
nentafterconvolution.Wealsoneedtoperformthisconvolutionintheproduct
spaceofthedomainandtherangeoftheinputsignal.ObservingEquations(1),
thenonlinearitycomesfromthedivisionbyW
bf
andfromthedependencyon
thepixelintensitiesthroughG
σ
r
(|I
p
?I
q
|).Westudyeachpointseparatelyand
isolatetheminthecomputationflow.
4
3.1HomogeneousIntensity
AdirectsolutiontohandlethedivisionistomultiplybothsidesofEquation(1a)
byW
bf
p
.Thetwoequationsarethenalmostsimilar.Weunderlinethispointby
rewritingEquations(1)usingtwo-dimensionalvectors:
?
?
W
bf
p
I
bf
p
W
bf
p
?
?
=
summationdisplay
q∈S
G
σ
s
(||p?q||)G
σ
r
(|I
p
?I
q
|)
?
?
I
q
1
?
?
(2)
Tomaintainthepropertythatthebilateralfilterisaweightedmean,weintro-
duceafunctionWwhosevalueis1everywhere:
?
?
W
bf
p
I
bf
p
W
bf
p
?
?
=
summationdisplay
q∈S
G
σ
s
(||p?q||)G
σ
r
(|I
p
?I
q
|)
?
?
W
q
I
q
W
q
?
?
(3)
Byassigningacouple(W
q
I
q
,W
q
)toeachpixelq,weexpressthefiltered
pixelsaslinearcombinationsoftheiradjacentpixels.Ofcourse,wehavenot
“removed”thedivisionsincetoaccesstheactualvalueoftheintensity,the
firstcoordinate(WI)hasstilltobedividedbythesecondone(W).Thiscanbe
comparedwithhomogeneouscoordinatesusedinprojectivegeometry.Addingan
extracoordinatetoourdatamakesmostofthecomputationpipelinecomputable
withlinearoperations;adivisionismadeonlyatthefinalstage.Inspiredbythis
parallel,wecallthecouple(WI,W)thehomogeneousintensity.
AlthoughEquation(3)isalinearcombination,thisdoesnotdefinealinear
filteryetsincetheweightsdependontheactualvaluesofthepixels.Thenext
sectionaddressesthisissue.
3.2TheBilateralFilterasaConvolution
IfweignorethetermG
σ
r
(|I
p
?I
q
|),Equation(3)isaclassicalconvolutionbya
Gaussiankernel:(W
bf
I
bf
,W
bf
)=G
σ
s
?(WI,W).Buttherangeweightdepends
onI
p
?I
q
andthereisnosummationonI.Toovercomethispoint,weintroduce
anadditionaldimensionζandsumoverit.WiththeKroneckersymbolδ(ζ)
(1ifζ=0,0otherwise)andRtheintervalonwhichtheintensityisdefined,we
rewriteEquation(3)using
bracketleftbig
δ(ζ?I
q
)=1
bracketrightbig
?
bracketleftbig
ζ=I
q
bracketrightbig
:
?
?
W
bf
p
I
bf
p
W
bf
p
?
?
=
summationdisplay
q∈S
summationdisplay
ζ∈R
G
σ
s
(||p?q||)G
σ
r
(|I
p
?ζ|)δ(ζ?I
q
)
?
?
W
q
I
q
W
q
?
?
(4)
Equation(4)isasumovertheproductspaceS×R.Wenowfocusonthisspace.
WeuselowercasenamesforthefunctionsdefinedonS×R.TheproductG
σ
s
G
σ
r
definesaGaussiankernelg
σ
s
,σ
r
onS×R:
g
σ
s
,σ
r
:(x∈S,ζ∈R)mapsto→G
σ
s
(||x||)G
σ
r
(|ζ|)(5)
FromtheremainingpartofEquation(4),webuildtwofunctionsiandw:
i:(x∈S,ζ∈R)mapsto→I
x
(6a)
w:(x∈S,ζ∈R)mapsto→δ(ζ?I
x
)W
x
(6b)
5
Thefollowingrelationsstemdirectlyfromthetwopreviousdefinitions:
I
x
=i(x,I
x
)andW
x
=w(x,I
x
)and?ζnegationslash=I
x
,w(x,ζ)=0(7)
ThenEquation(4)isrewrittenas:
?
?
W
bf
p
I
bf
p
W
bf
p
?
?
=
summationdisplay
(q,ζ)∈S×R
g
σ
s
,σ
r
(p?q,I
p
?ζ)
?
?
w(q,ζ)i(q,ζ)
w(q,ζ)
?
?
(8)
Theaboveformulacorrespondstothevalueatpoint(p,I
p
)ofaconvolution
betweeng
σ
s
,σ
r
andthetwo-dimensionalfunction(wi,w):
?
?
W
bf
p
I
bf
p
W
bf
p
?
?
=
?
?
g
σ
s
,σ
r
?
?
?
wi
w
?
?
?
?
(p,I
p
)(9)
Accordingtotheaboveequation,weintroducethefunctionsi
bf
andw
bf
:
(w
bf
i
bf
,w
bf
)=g
σ
s
,σ
r
?(wi,w)(10)
Thus,wehavereachedourgoal.Thebilateralfilterisexpressedasaconvolution
followedbynonlinearoperations:
linear:(w
bf
i
bf
,w
bf
)=g
σ
s
,σ
r
?(wi,w)(11a)
nonlinear:I
bf
p
=
w
bf
(p,I
p
)i
bf
(p,I
p
)
w
bf
(p,I
p
)
(11b)
Thenonlinearsectionisactuallycomposedoftwooperations.Thefunctions
w
bf
i
bf
andw
bf
areevaluatedatpoint(p,I
p
).Wenamethisoperationslicing.
Thesecondnonlinearoperationisthedivision.Inourcase,slicinganddivision
commutei.e.theresultisindependentoftheirorderbecauseg
σ
s
,σ
r
ispositive
andwvaluesare0and1,whichensuresthatw
bf
ispositive.
3.3Intuition
Togainmoreintuitionaboutourformulationofthebilateralfilter,wepropose
aninformaldescriptionoftheprocessbeforediscussingfurtheritsconsequences.
ThespatialdomainSisaclassicalxyimageplaneandtherangedomain
Risasimpleaxislabelledζ.Thewfunctioncanbeinterpretedasthe“plot
inthexyζspaceofζ=I(x,y)”i.e.wisnulleverywhereexceptonthepoints
(x,y,I(x,y))whereitisequalto1.Thewiproductissimilartow.Insteadof
usingbinaryvalues0or1to“plotI”,weuse0orI(x,y)i.e.itisaplotwitha
penwhosebrightnessequalstheplottedvalue.
Thenusingthesetwofunctionswiandw,thebilateralfilteriscomputedas
follows.First,we“blur”wiandwi.e.weconvolvewiandwwithaGaussian
definedonxyζ.Thisresultsinthefunctionsw
bf
i
bf
andw
bf
.Foreachpoint
6
x
ζ
x
ζ
x
ζ
x
ζ
x
ζ
samplinginthexζspace
space(x)
range(
ζ
)
Gaussianconvolution
division
slicing
0
0.2
0.4
0.6
0.8
1
020406080100120
w
w
bf
i
bf
w
bf
wi
0
0.2
0.4
0.6
0.8
1
020406080100120
Fig.1.Ourcomputationpipelineappliedtoa1Dsignal.Theoriginaldata(toprow)
arerepresentedbyatwo-dimensionalfunction(wi,w)(secondrow).Thisfunctionis
convolvedwithaGaussiankerneltoform(w
bf
i
bf
,w
bf
)(thirdrow).Thefirstcomponent
isthendividedbythesecond(fourthrow,blueareaisundefinedbecauseofnumerical
limitation,w
bf
≈0).Thenthefinalresult(lastrow)isextractedbysamplingthe
formerresultatthelocationoftheoriginaldata(showninredonthefourthrow).
ofthexyζspace,wecomputei
bf
(x,y,ζ)bydividingw
bf
(x,y,ζ)i
bf
(x,y,ζ)by
w
bf
(x,y,ζ).Thefinalstepistogetthevalueofthepixel(x,y)ofthefiltered
imageI
bf
.Thisdirectlycorrespondstothevalueofi
bf
at(x,y,I(x,y))which
isthepointwheretheinputimageIwas“plotted”.Figure1illustratesthis
processonasimple1Dimage.
4FastApproximation
WehaveshownthatthebilateralfiltercanbeinterpretedasaGaussianfilterin
aproductspace.Ouraccelerationschemedirectlyfollowsfromthefactthatthis
operationisalow-passfilter.(w
bf
i
bf
,w
bf
)isthereforeessentiallyaband-limited
functionwhichiswellapproximatedbyitslowfrequencies.
Usingthesamplingtheorem[20](p.35),itissu?cienttosamplewitharate
atleasttwiceshorterthanthesmallestwavelengthconsidered.Inpractice,we
downsample(wi,w),performtheconvolution,andupsampletheresult:
(w
↓
i
↓
,w
↓
)=downsample(wi,w)(12a)
(w
bf
↓
i
bf
↓
,w
bf
↓
)=g
σ
s
,σ
r
?(w
↓
i
↓
,w
↓
)(12b)
(w
bf
↓↑
i
bf
↓↑
,w
bf
↓↑
)=upsample(w
bf
↓
i
bf
↓
,w
bf
↓
)(12c)
7
Therestofthecomputationremainsthesameexceptthatwesliceanddivide
(w
bf
↓↑
i
bf
↓↑
,w
bf
↓↑
)insteadof(w
bf
i
bf
,w
bf
),usingthesame(p,I
p
)points.Sinceslicing
occursatpointswherew=1,itguaranteesw
bf
≥g
σ
s
,σ
r
(0),whichensuresthat
wedonotdividebysmallnumbersthatwoulddegradeourapproximation.
Weusebox-filteringfortheprefilterofthedownsampling(a.k.a.average
downsampling),andlinearupsampling.Whilethesefiltersdonothaveperfect
frequencyresponses,theyo?ermuchbetterperformancesthanschemessuchas
tri-cubicfilters.
4.1Evaluation
Toevaluatetheerrorinducedbyourapproximation,wecomparetheresult
I
bf
↓↑
fromthefastalgorithmtoI
bf
obtainedfromEquations(1).Forthispur-
pose,wecomputethepeaksignal-to-noiseratio(PSNR)consideringR=[0;1]:
PSNR(I
bf
↓↑
)=?10log
10
parenleftbigg
1
|S|
summationtext
p∈S
vextendsingle
vextendsingle
vextendsingleI
bf
↓↑
(p)?I
bf
(p)
vextendsingle
vextendsingle
vextendsingle
2
parenrightbigg
.
Wehavechosenthreeimagesasdi?erentaspossibletocoverabroadspec-
trumofcontent.Weuse(seeFigure4):
–Anartificialimagewithvariousedgesandfrequencies,andwhitenoise.
–Anarchitecturalpicturestructuredalongtwomaindirections.
–Andaphotographofanaturalscenewithamorestochasticstructure.
Theboxdownsamplingandlinearupsamplingschemesyieldverysatisfying
resultswhilebeingcomputationallye?cient.Weexperimentedseveralsampling
rates(s
s
,s
r
)forS×R.Themeaningfulquantitiestoconsideraretheratios
parenleftBig
s
s
σ
s
,
s
r
σ
r
parenrightBig
thatindicatehowmanyhighfrequenciesweignorecomparedtothe
bandwidthofthefilterweapply.Smallratioscorrespondtolimitedapproxima-
tionsandhighratiostomoreaggressivedownsamplings.Aconsistentapproxi-
mationisasamplingrateproportionaltotheGaussianbandwidth(i.e.
s
s
σ
s
≈
s
r
σ
r
)
toachievesimilaraccuracyonthewholeS×Rdomain.Theresultsplottedin
Figure2showthatthisremarkisgloballyvalidinpractice.Acloserlookat
theplotsrevealsthatScanbeslightlymoreaggressivelydownsampledthanR.
Thisisprobablyduetothenonlinearitiesandtheanisotropyofthesignal.
0.40.20.10.050.025
64
32
16
8
4
0.40.20.10.050.025
64
32
16
8
4
0.40.20.10.050.025
64
32
16
8
4
intensitysampling
spacesampling
intensitysamplingintensitysampling
spacesamplingspacesampling
(a)artificial(b)architectural(c)natural
50dB
50dB
50dB
40dB
40dB
40dB
30dB30dB
30dB
Fig.2.Accuracyevaluation.Alltheimagesarefilteredwith(σ
s
=16,σ
r
=0.1).The
PSNRindBisevaluatedatvarioussamplingofSandR(greaterisbetter).Ourap-
proximationschemeismorerobusttospacedownsamplingthanrangedownsampling.
Itisalsoslightlymoreaccurateonstructuredscenes(a,b)thanstochasticones(c).
8
15
20
25
30
35
40
45
50
55
60
110100
PSNR(indB)
Durand-Dorseyapproximation
ourapproximation
0.4
time(ins)
0.40.20.10.050.025
64
32
16
8
4
intensitysampling
spacesampling
40s
20s
10s
5s
2s
1s
0.8s
0.6s
Fig.3.Left:Runningtimesonthearchitecturalpicturewith(σ
s
=16,σ
r
=0.1).
ThePSNRisolinesareplottedingray.Exactcomputationtakesabout1h.?Right:
Accuracy-versus-timecomparison.Bothmethodsaretestedonthearchitecturalpicture
(1600×1200)withthesamesamplingratesofS×R(fromlefttoright):(4;0.025)(8;0.05)
(16;0.1)(32,0.2)(64,0.4).
Figure3-leftshowstherunningtimesforthearchitecturalpicturewiththe
samesettings.Intheory,thegainfromspacedownsamplingshouldbetwicethe
onefromrangedownsamplingsinceSistwo-dimensionalandRone-dimensional.
Inpractice,thenonlinearitiesandcachingissuesinduceminordeviations.Com-
biningthisplotwiththePSNRplot(ingrayundertherunningtimes)allows
forselectingthebestsamplingparametersforagivenerrortoleranceoragiven
timebudget.Asasimpleguideline,usingsamplingstepsequaltoσ
s
andσ
r
pro-
duceresultswithoutvisualdi?erencewiththeexactcomputation(seeFig.4).
Ourschemeachievesaspeed-upoftwoordersofmagnitude:Directcomputation
ofEquations(1)lastsaboutonehourwhereasourapproximationrequiresone
second.Thisdramaticimprovementopensavenuesforinteractiveapplications.
4.2ComparisonwiththeDurand-Dorseyspeed-up
DurandandDorseyalsodescribealinearapproximationofthebilateralfilter[3].
UsingevenlyspacedintensityvaluesI
1
..I
n
thatcoverR,theirschemecanbe
summarizedas(forconvenience,wealsonameG
σ
s
the2DGaussiankernel):
ι
↓
=downsample(I)[imagedownsampling](13a)
?ζ∈{1..n}ω
↓ζ
(p)=G
σ
r
(|ι
↓
(p)?I
ζ
|)[rangeweightevaluation](13b)
?ζ∈{1..n}ωι
↓ζ
(p)=ω
↓ζ
(p)ι
↓
(p)[intensitymultiplication](13c)
?ζ∈{1..n}(ωι
bf
↓ζ
,ω
bf
↓ζ
)=G
σ
s
?(ωι
↓ζ
,ω
↓ζ
)[convolutiononS](13d)
?ζ∈{1..n}ι
bf
↓ζ
=ωι
bf
↓ζ
/ω
bf
↓ζ
[normalization](13e)
?ζ∈{1..n}ι
bf
↓↑ζ
=upsample(ι
bf
↓ζ
)[layerupsampling](13f)
I
DD
(p)=interpolation(ι
bf
↓↑ζ
)(p)[nearestlayerinterpolation](13g)
Withoutdownsampling(i.e.{I
i
}=RandSteps13a,fignored),theDurand-
DorseyschemeisequivalenttooursbecauseSteps13b,c,dcorrespondtoacon-
volutiononS×R.Indeed,Step(13b)computesthevaluesofaGaussiankernel
onR.Step(13c)actuallyevaluatestheconvolutiononR,consideringthatι
↓
(p)
istheonlynonzerovalueontheζaxis.WithStep(13d),theconvolutiononS,
thesethreestepsperforma3DconvolutionusingaseparationbetweenRandS.
9
Thedi?erencescomesfromthedownsamplingapproach.DurandandDorsey
interleavelinearandnonlinearoperations:Thedivisionisdoneaftertheconvo-
lution13dbutbeforetheupsampling13f.Thereisnosimpletheoreticalbaseto
estimatetheerror.Moreimportantly,theDurand-Dorseystrategyissuchthat
theintensityιandtheweightωarefunctionsdefinedonSonly.Agivenpixel
hasonlyoneintensityandoneweight.Afterdownsampling,bothsidesofthe
discontinuitymayberepresentedbythesamevaluesofιandω.Thisisapoor
representationofthediscontinuitiessincetheyinherentlyinvolveseveralvalues.
Incomparison,wedefinefunctionsonS×R.ForagivenimagepointinS,we
canhandleseveralvaluesontheRdomain.TheadvantageofworkinginS×R
isthatthischaracteristicisnotalteredbydownsampling.Itisthemajorreason
whyourschemeismoreaccuratethantheDurand-Dorseytechnique,especially
ondiscontinuities.
Figure3-rightshowstheprecisionachievedbybothapproachesrelativelyto
theirrunningtime,andFigure4illustratestheirvisualdi?erences.Thereisno
gaininextremedownsamplingsincenonlinearitiesarenomorenegligible.Both
approachesalsohaveaplateauintheiraccuracyi.e.beyondacertainpoint,pre-
cisiongainsincreaseslowlywithsamplingrefinementbutoursreachesahigher
accuracy(≈55dBcomparedto≈40dB).Inaddition,forthesamerunningtime,
ourapproachisalwaysmoreaccurate(exceptforextremedownsampling).
4.3Implementation
AlltheexperimentshavebeendoneanIntelXeon2.8GHzusingthesamecode
baseinC++.WehaveimplementedtheDurand-Dorseytechniquewiththe
samelibrariesasourtechnique.2Dand3DconvolutionsaremadeusingFFT.
Thedomainsarepaddedwithzerosover2σtoavoidcross-boundaryartefacts.
Therearenoothersignificantoptimizationstoavoidbiasinthecomparisons.
Aproductionimplementationcouldthereforebeimprovedwithtechniquessuch
asseparableconvolution.Ourcodeispubliclyavailableonourwebpage.The
softwareisopen-source,undertheMITlicense.
5Discussion
DimensionalityOurseparationintolinearandnonlinearpartscomesatthe
costoftheadditionalζdimension.Onehastobecarefulbeforeincreasingthe
dimensionalityofaproblemsincetheincurredperformanceoverheadmayexceed
thegains,restrictingourstudytoatheoreticaldiscussion.Wehavehowever
demonstratedthatthisformalismallowsforacomputationschemethatisseveral
ordersofmagnitudefasterthanastraightforwardapplicationofEquation(1).
ThisadvocatesperformingthecomputationintheS×Rspaceinsteadofthe
imageplane.Inthisrespect,ourapproachcanbecomparedtothelevelsets[21]
whichalsodescribeabettercomputationschemeusingahigherdimensionspace.
Notethatusingthetwo-dimensionalhomogeneousintensitydoesnotincrease
thedimensionalitysinceEquation(1)alsocomputestwofunctions:W
bf
andI
bf
.
ComparisonwithGeneralizedIntensityBarashalsousespointsintheS×R
spacethathenamesgeneralizedintensities[13].Ourtwoapproacheshavein
commontheglobalmeaningofSandR:Theformerisrelatedtopixelpositions
10
sampling(σ
s
,σ
r
)(4,0.025)(8,0.05)(16,0.1)(32,0.2)(64,0.4)
downsampling1.3s0.23s0.09s0.07s0.06s
convolution63s2.8s0.38s0.02s0.01s
nonlinearity0.48s0.47s0.46s0.47s0.46s
Table1.Timeusedbyeachstepatdi?erentsamplingratesofthearchitecturalimage.
Upsamplingisincludedinthenonlinearitytimebecauseourimplementationcomputes
i
bf
↓↑
onlyatthe(x,I
x
)pointsratherthanupsamplingthewholeS×Rspace.
andthelattertopixelintensities.Itisneverthelessimportanttohighlightthe
di?erences.BarashhandlesS×Rtocomputedistances.Thus,hecanexpress
thedi?erencebetweenadaptivesmoothingandbilateralfilteringasadi?erence
ofdistancedefinitions.Buttheactualcomputationremainsthesame.Ouruse
ofS×Rismoreinvolved.Wenotonlymanipulatepointsinthisspacebut
alsodefinefunctionsandperformconvolutionsandslicing.Anotherdi?erence
isourdefinitionoftheintensitythroughafunction(iori
bf
).Barashassociates
directlytheintensityofapointtoitsRcomponentwhereasinourframework,
theintensityofapoint(x,ζ)isnotitsζcoordinatee.g.ingenerali
bf
(x,ζ)negationslash=ζ.
Inaddition,ourapproachleadstoamoree?cientimplementation.
ComplexityOneoftheadvantageofourseparationisthattheconvolutionis
themostcomplexpartofthealgorithm.Using|·|forthecardinalofaset,the
convolutioncanbedoneinO(|S||R|log(|S||R|))withfastFouriertransform
andmultiplicationinthefrequencydomain.Thentheslicingandthedivision
aredonepixelbypixel.Thustheyarelinearintheimagesizei.e.O(|S|).Hence,
thealgorithmcomplexityisdominatedbytheconvolution.Thisresultisverified
inpracticeasshownbyTable1.Theconvolutiontimerapidlyincreasesasthe
samplingbecomesfiner.Thisvalidatesourchoiceoffocussingontheconvolution.
6Conclusions
Wehavepresentedafastapproximationtechniqueofthebilateralfilterbased
onasignalprocessinginterpretation.Fromatheoreticalpointofview,wehave
introducedthenotionofhomogeneousintensityanddemonstratedanewap-
proachofthespace-intensitydomain.Webelievethattheseconceptscanbe
appliedbeyondbilateralfiltering,andwehopethatthesecontributionswillin-
spirenewstudies.Fromapracticalpointofview,ourapproximationtechnique
yieldsresultsvisuallysimilartotheexactcomputationwithinteractiverun-
ningtimes.Thistechniquepavesthewayforinteractiveapplicationsrelyingon
qualityimagesmoothing.
FutureWorkOurstudytranslatesalmostdirectlytohigherdimensionaldata
(e.g.colorimagesorvideos).Analyzingtheperformanceinthesecaseswill
providevaluablestatistics.ExploringdeeperthefrequencystructureoftheS×R
domainseemsanexcitingresearchdirection.
AcknowledgementWethankSoonminBae,SamuelHornus,ThouisJones,and
SaraSufortheirhelpwiththepaper.ThisworkwassupportedbyaNational
ScienceFoundationCAREERaward0447561“TransientSignalProcessingfor
11
RealisticImagery,”anNSFGrantNo.0429739“ParametricAnalysisandTrans-
ferofPictorialStyle,”agrantfromRoyalDutch/ShellGroupandtheOxygen
consortium.Fr′edoDurandacknowledgesaMicrosoftResearchNewFacultyFel-
lowship,andSylvainPariswaspartiallysupportedbyaLavoisierFellowship
fromtheFrench“Minist`eredesA?aires
′
Etrang`eres.”
References
1.Tomasi,C.,Manduchi,R.:Bilateralfilteringforgrayandcolorimages.In:Proc.
ofInternationalConferenceonComputerVision,IEEE(1998)839–846
2.Oh,B.M.,Chen,M.,Dorsey,J.,Durand,F.:Image-basedmodelingandphoto
editing.In:Proc.ofSIGGRAPHconference,ACM(2001)
3.Durand,F.,Dorsey,J.:Fastbilateralfilteringforthedisplayofhigh-dynamic-range
images.ACMTrans.onGraphics21(2002)Proc.ofSIGGRAPHconference.
4.Eisemann,E.,Durand,F.:Flashphotographyenhancementviaintrinsicrelighting.
ACMTrans.onGraphics23(2004)Proc.ofSIGGRAPHconference.
5.Petschnigg,G.,Agrawala,M.,Hoppe,H.,Szeliski,R.,Cohen,M.,Toyama,K.:
Digitalphotographywithflashandno-flashimagepairs.ACMTrans.onGraphics
23(2004)Proc.ofSIGGRAPHconference.
6.Jones,T.R.,Durand,F.,Desbrun,M.:Non-iterative,feature-preservingmesh
smoothing.ACMTrans.onGraphics22(2003)Proc.ofSIGGRAPHconference.
7.Fleishman,S.,Drori,I.,Cohen-Or,D.:Bilateralmeshdenoising.ACMTrans.on
Graphics22(2003)Proc.ofSIGGRAPHconference.
8.Wong,W.C.K.,Chung,A.C.S.,Yu,S.C.H.:Trilateralfilteringforbiomedicalim-
ages.In:Proc.ofInternationalSymposiumonBiomedicalImaging,IEEE(2004)
9.Bennett,E.P.,McMillan,L.:Videoenhancementusingper-pixelvirtualexposures.
ACMTrans.onGraphics24(2005)845–852Proc.ofSIGGRAPHconference.
10.Paris,S.,Brice?no,H.,Sillion,F.:Captureofhairgeometryfrommultipleimages.
ACMTrans.onGraphics23(2004)Proc.ofSIGGRAPHconference.
11.Elad,M.:Onthebilateralfilterandwaystoimproveit.IEEETrans.OnImage
Processing11(2002)1141–1151
12.Smith,S.M.,Brady,J.M.:SUSAN–anewapproachtolowlevelimageprocessing.
InternationalJournalofComputerVision23(1997)45–78
13.Barash,D.:Afundamentalrelationshipbetweenbilateralfiltering,adaptive
smoothingandthenonlineardi?usionequation.IEEETrans.onPatternAnalysis
andMachineIntelligence24(2002)844
14.Buades,A.,Coll,B.,Morel,J.M.:NeighborhoodfiltersandPDE’s.Technical
Report2005-04,CMLA(2005)
15.Yaroslavsky,L.P.:DigitalPictureProcessing.SpringerVerlag(1985)
16.Huber,P.J.:RobustStatistics.Wiley-Interscience(1981)
17.Hampel,F.R.,Ronchetti,E.M.,Rousseeuw,P.M.,Stahel,W.A.:RobustStatistics
–TheApproachBasedonInfluenceFunctions.WileyInterscience(1986)
18.Mr′azek,P.,Weickert,J.,Bruhn,A.:OnRobustEstimationandSmoothing
withSpatialandTonalKernels.In:GeometricPropertiesfromIncompleteData.
Springer(toappear)
19.Choudhury,P.,Tumblin,J.E.:Thetrilateralfilterforhighcontrastimagesand
meshes.In:Proc.ofEurographicsSymposiumonRendering.(2003)
20.Smith,S.:DigitalSignalProcessing.Newnes(2002)
21.Osher,S.,Sethian,J.A.:Frontspropagatingwithcurvature-dependentspeed:Al-
gorithmsbasedonHamilton-Jacobiformulations.J.ofComp.Physics.(1988)
12
originalimage
exactbilateralfilter
approximatedbilateralfilter
dif
ference
Durand-Dorsey
approximation
artificialarchitecturalnatural
Fig.4.Wehavetestedourapproximatedschemeonthreeimages(firstrow):anar-
tificialimage(512×512)withdi?erenttypesofedgesandawhitenoiseregion,an
architecturalpicture(1600×1200)withstrongandorientedfeatures,andanatural
photograph(800×600)withmorestochastictextures.Forclarity,wepresentrepresen-
tativeclose-ups(secondrow).Fullresolutionimagesareavailableonourwebsite.Our
approximationproducesresults(fourthrow)visuallysimilartotheexactcomputation
(thirdrow).Acolorcodedsubtraction(fifthrow)revealssubtledi?erencesattheedges
(red:negative,black:0,andblue:positive).Incomparison,theDurand-Dorseyapprox-
imationintroduceslargevisualdiscrepancies:thedetailsarewashedout(bottomrow).
Allthefiltersarecomputedwithσ
s
=16andσ
r
=0.1.Ourfilterusesasamplingrate
of(16,0.1).ThesamplingrateoftheDurand-Dorseyfilterischoseninordertoachieve
thesame(orslightlysuperior)runningtime.Thus,thecomparisonisdonefairly,using
thesametimebudget.
|
|
