Felzenszwalb Computer Science Department University of Chicago pffcsuchicagoedu Daniel P Huttenlocher Computer Science Department Cornell University dphcscornelledu Abstract Markov random 64257eld models provide a robust and uni64257ed framew ID: 8034
Download Pdf The PPT/PDF document "Efcient Belief Propagation for Early Vis..." is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.
EfcientBeliefPropagationforEarlyVisionPedroF.FelzenszwalbComputerScienceDepartment,UniversityofChicagopff@cs.uchicago.eduDanielP.HuttenlocherComputerScienceDepartment,CornellUniversitydph@cs.cornell.edu1 AbstractMarkovrandomeldmodelsprovidearobustanduniedframeworkforearlyvisionproblemssuchasstereoandimagerestoration.Inferencealgorithmsbasedongraphcutsandbeliefpropa-gationhavebeenfoundtoyieldaccurateresults,butdespiterecentadvancesareoftentooslowforpracticaluse.Inthispaperwepresentsomealgorithmictechniquesthatsubstantiallyimprovetherunningtimeoftheloopybeliefpropagationapproach.Oneofthetechniquesreducesthecomplex-ityoftheinferencealgorithmtobelinearratherthanquadraticinthenumberofpossiblelabelsforeachpixel,whichisimportantforproblemssuchasimagerestorationthathavealargelabelset.Anothertechniquespeedsupandreducesthememoryrequirementsofbeliefpropagationongridgraphs.Athirdtechniqueisamulti-gridmethodthatmakesitpossibletoobtaingoodresultswithasmallxednumberofmessagepassingiterations,independentofthesizeoftheinputimages.Takentogetherthesetechniquesspeedupthestandardalgorithmbyseveralordersofmagnitude.Inpracticeweobtainresultsthatareasaccurateasthoseofotherglobalmethods(e.g.,usingtheMiddleburystereobenchmark)whilebeingnearlyasfastaspurelylocalmethods.1IntroductionOverthepastfewyearstherehavebeenexcitingadvancesinthedevelopmentofalgorithmsforsolvingearlyvisionproblemssuchasstereoandimagerestorationusingMarkovrandomeld(MRF)models.WhiletheMRFframeworkyieldsanoptimizationproblemthatisNPhard,goodapproximationtechniquesbasedongraphcuts[4]andonbeliefpropagation[12,10]havebeendevelopedanddemonstratedforproblemssuchasstereoandimagerestoration.Thesemethodsaregoodbothinthesensethatthelocalminimatheyndareminimaoverlargeneighborhoods,andinthesensethattheyproducehighlyaccurateresultsinpractice.Acomparisonbetweenthetwodifferentapproachesforthecaseofstereoisdescribedin[11].Despitethesesubstantialadvanceshowever,boththegraphcutsandbeliefpropagationap-proachesarestillcomputationallydemandingwhencomparedtolocalmethodsthatarenotat-temptingtosolveMRF-basedformulations.ThusoneisfacedwithchoosingbetweentheMRF-2 basedmethods,whichproducegoodresultsbutarerelativelyslow,andlocalmethodswhichpro-ducesubstantiallypoorerresultsbutarefast.Inthispaperwepresentseveralalgorithmictech-niquesthatsubstantiallyimprovetherunningtimeofbeliefpropagation(BP)forsolvingearlyvisionproblems.Takentogetherthesetechniquesspeedupthestandardalgorithmbyseveralor-dersofmagnitude,makingitsrunningtimecompetitivewithlocalmethods.InthecaseofstereoweobtainresultswithacomparabledegreeofaccuracytostandardBPorgraphcutsalgorithmsinlessthanonesecondperimagepair.Thedifferencesareevenmorepronouncedforproblemssuchasimagerestorationwheretherearearelativelylargenumberofpossiblelabelsforeachpixel.Thegeneralframeworkfortheproblemsweconsidercanbedenedasfollows(wefollowthenotationin[4]).LetPbethesetofpixelsinanimageandLbeanitesetoflabels.Thelabelscorrespondtoquantitiesthatwewanttoestimateateachpixel,suchasdisparitiesorintensities.Alabelingfassignsalabelfp2Ltoeachpixelp2P.Weassumethatthelabelsshouldvaryslowlyalmosteverywherebutmaychangedramaticallyatsomeplacessuchaspixelsalongobjectboundaries.Thequalityofalabelingisgivenbyanenergyfunction,E(f)=Xp2Pp(fp)+X(p;q)2NW(fp;fq);whereNarethe(undirected)edgesinthefour-connectedimagegridgraph.p(fp)isthecostofassigninglabelfptopixelp,andisreferredtoasthedatacost.W(fp;fq)measuresthecostofassigninglabelsfpandfqtotwoneighboringpixels,andisnormallyreferredtoasthediscontinuitycost.Findingalabelingthatminimizesthisenergycorrespondstothemaximumaposteriori(MAP)estimationproblemforanappropriatelydenedMRF(see[2,8]).Inlow-levelcomputervisionproblemsthediscontinuitycostWisgenerallybasedonthedifferencebetweenlabels,ratherthanontheiractualvalues.Forexample,instereoandimagerestorationthelabelscorrespondtothepossibledisparitiesorintensityvaluesrespectively,andthecostofassigningapairoflabelstoneighboringpixelsisbasedonthedegreeofdifferencebetweenthelabels.ThusinthispaperweconsiderthecasewhereW(fp;fq)=V(fp fq),yieldinganenergyminimizationproblemoftheform,E(f)=Xp2Pp(fp)+X(p;q)2NV(fp fq):(1)3 2LoopyBeliefPropagation:Max-ProductWestartbybrieyreviewingtheBPapproachforperforminginferenceonMarkovrandomelds(e.g.,see[12]).Firstweconsiderthemax-productalgorithm,whichcanbeusedtoapproximatetheMAPsolutiontoMRFproblems.Normallythistechniqueisdenedintermsofprobabilitydistri-butions,butanequivalentcomputationcanbeperformedwithnegativelogprobabilities,wherethemax-productbecomesamin-sum.Weusethisformulationbecauseitislesssensitivetonumericalartifacts,andbecauseitdirectlycorrespondstotheenergyfunctiondenitioninequation(1).Themax-productBPalgorithmworksbypassingmessagesaroundthegraphdenedbythefour-connectedimagegrid.Themethodisiterative,withmessagesfromallnodesbeingpassedinparallel.Eachmessageisavectorofdimensiongivenbythenumberofpossiblelabels,k.Letmtp!qbethemessagethatnodepsendstoaneighboringnodeqatiterationt.Whenusingnegativelogprobabilitiesallentriesinm0p!qareinitializedtozero,andateachiterationnewmessagesarecomputedinthefollowingway,mtp!q(fq)=minp0@V(fp fq)+p(fp)+Xs2N(p)nqmt 1s!p(fp)1A(2)whereN(p)nqdenotestheneighborsofpotherthanq.AfterTiterationsabeliefvectoriscom-putedforeachnode,q(fq)=q(fq)+Xp2N(q)mTp!q(fq):Finally,thelabelfqthatminimizesq(fq)individuallyateachnodeisselected.Thestandardimplementationofthismessagepassingalgorithmonthegridgraphrunsin(nk2T)time,wherenisthenumberofpixelsintheimage,kisthenumberofpossiblelabelsforeachpixelandTisthenumberofiterations.Ittakes(k2)timetocomputeeachmessageandthereare(n)messagestobecomputedineachiteration.Weconsiderthreedifferenttechniquesforspeedingupthisstandardbeliefpropagationalgo-rithm.Firstwenotethateachmessageupdatecanbeexpressedasaminconvolution,andmoreoverthatwiththediscontinuitycostscommonlyusedinearlyvisionthisminconvolutioncanbecom-putedin(k)time.Oursecondresultshowsthatforthegridgraph(andanybipartitegraph)4 essentiallythesamebeliefsasthosedenedabovecanbeobtainedusingonlyhalfasmanymes-sageupdates.Besidesyieldingaspeedupthistechniquealsomakesitpossibletocomputethemessagesinplace,usinghalfasmuchmemoryasthenormalalgorithm.Thisisimportantbe-causeBPhashighmemoryrequirements,storingmultipledistributionsateachpixel.Finallywepresentamulti-gridapproachforperformingBPinacoarse-to-nemanner.Inthisapproachthenumberofmessagepassingiterations,T,canbeasmallconstant,becauselongrangeinteractionsarecapturedbyshortpathsincoarsegridgraphs.Combiningthethreetechniquestogetheryieldsanalgorithmthatrunsin(nk)timeoverallandisveryfastinpractice.Incontrastthestandardalgorithmsummarizedaboverequires(nk2)timeperiterationandthenumberofiterationsTisgenerally(n1=2)toallowforinformationtopropagateacrosstheimage.Ourexperimentalresultsareasaccurateasthoseobtainedwhenusingstandardmax-productBPorgraphcutsalgorithmstominimizeenergyfunctionsoftheforminequation(1).Inthecaseofstereowequantifythisusingthebenchmarkin[9].3ComputingaMessageUpdateinLinearTimeThissectioncoverstherstofthethreetechniques,whichreducesthetimerequiredtocomputeasinglemessageupdatefrom(k2)to(k)formostlow-levelvisionapplications.Wecanre-writeequation(2)as,mtp!q(fq)=minp(V(fp fq)+h(fp));(3)whereh(fp)=p(fp)+Pmt 1s!p(fp).Thestandardwayofcomputingthismessageistoexplicitlyminimizeoverfpforeachchoiceoffq,whichtakes(k2)timewherekisthenumberoflabels.Theformofequation(3)iscommonlyreferredtoasaminconvolution.Thisoperationisanalogoustothestandarddiscreteconvolutionoperation,howeverinastandardconvolutionthesumwouldbereplacedbyaproductandtheminwouldbereplacedbyasum.WhilestandarddiscreteconvolutionscanbeefcientlycomputedusingtheFFT,nosuchgeneralresultisknownforminconvolutions.However,forthecostfunctionsV(fp fq)commonlyusedincomputervisionweshowthattheminconvolution,andthustheBPmessageupdates,canoftenbecomputed5 in(k)time.Suchlineartimemethodsareparticularlyimportantforproblemssuchasimagerestorationwherethenumberoflabels,k,canbeinthehundredsormore.Notethatthesemethodsdonotmakeanyapproximations;theycomputeexactlythesameresultasthe(k2)bruteforcealgorithm.3.1PottsModelWestartbyconsideringasimplemeasureofthedifferencebetweenlabels,thePottsmodel[4],whichcapturestheassumptionthatlabelingsshouldbepiecewiseconstant.Thismodelconsidersonlytheequalityorinequalityoflabels.Forequallabelsthecostiszero,whilefordifferentlabelsthecostisapositiveconstant,V(x)=8:0ifx=0dotherwiseWiththiscostfunctiontheminconvolutioninequation(3)canbeexpressedinaformwheretheminimizationoverfpcanbeperformedonce,independentofthevalueoffq,mtp!q(fq)=minh(fq);minph(fp)+d:Separatingtheminimizationoverfpinthismannerreducesthetimenecessarytocomputeames-sageto(k).Firstwecomputeminph(fp),andthenusethattocomputethemessagevalueforeachfqinconstanttime.Notethatthisideastillapplieswheninsteadofasingleconstantdthereisaconstantdpqforeachedgeinthegraph.Thisisusefulwhentheresultofsomeotherprocess,suchasedgedetectionorsegmentation,suggeststhatdiscontinuitiesshouldbepenalizedmoreorlessfordifferentpairsofpixels.3.2LinearModelNowweconsiderthecasewherethecostfunctionVisbasedonthemagnitudeofthedifferencebetweenlabelsfpandfq.Onecommonsuchfunctionisthetruncatedlinearmodel,wherethecostincreaseslinearlybasedonthedistancebetweenthelabelsfpandfquptosomelevel.Inordertoallowforlargediscontinuitiesinthelabelingthecostfunctionstopsgrowingafterthedifference6 becomeslarge.Forinstance,V(x)=min(cx;d);(4)wherecistherateofincreaseinthecost,anddcontrolswhenthecoststopsincreasing.AsimilarcostfunctionwasusedinaBPapproachtostereo[10],althoughratherthantruncatingthelinearcosttheyhaveafunctionthatchangessmoothlyfrombeingalmostlinearneartheorigintoaconstantvalueasthecostincreases.Werstconsiderthesimplerproblemofapurelinearcostwithouttruncation.Substitutingintoequation(3)yields,mtp!q(fq)=minp(cfp fq+h(fp)):(5)Onecanenvisionthelabelsasbeingembeddedinagrid.Notethatthisisagridoflabelsandisnotrelatedtotheimagegrid.Forinstanceitisaone-dimensionalgridofdisparitylabelsinthecaseofstereoandaone-dimensionalgridofintensitylabelsinthecaseofimagerestoration.Theminimizationin(5)canthenbeseenasthelowerenvelopeofkupwardfacingconesofslopecrootedat(fp;h(fp))foreachgridpointfp.Theone-dimensionalcaseisillustratedinFigure1.Thislowerenvelopecalculationissimilartothatperformedincomputingadistancetransform(e.g.,[3]).Forthedistancetransformtheconesareplacedatheight0andoccuronlyatselectedvaluesratherthaneverygridpoint.Despitethesedifferences,thestandarddistancetransformalgorithmfrom[3]canbemodiedtocomputetheminconvolutionwithalinearcost.Itisstraightforwardtoverifythatthefollowingsimpletwo-passalgorithmcorrectlycomputesthemessageinequation(5)forthecasewherethelabelscorrespondtointegersf0;:::;k 1g.Firstweinitializethemessagevectorwiththevaluesofh,andthenupdateitsentriessequentially.Thisisdoneinplacesothatupdatesaffectoneanother,forfqfrom1tok 1:m(fq) min(m(fq);m(fq 1)+c):Thebackwardpassisanalogous,forfqfromk 2to0:m(fq) min(m(fq);m(fq+1)+c):7 Figure1:Anillustrationofthelowerenvelopeoffourconesinthecaseofone-dimensionallabels(e.g.stereodisparityorimagerestoration).Eachconeisrootedatlocation(fp;h(fp)).Thedarkerlineindicatesthelowerenvelope.ConsidertheexampleinFigure1.Theinitialvalueofmis(3;1;4;2).Withc=1,theforwardpassyields(3;1;2;2),andthebackwardpassyields(2;1;2;2).Thekeypropertythatallowsustousethisalgorithmisthatthelabelsareembeddedinagrid,andthediscontinuitycostisalinearfunctionofdistanceinthegrid.Ifthelabelsareembeddedinahigherdimensionalgrid(e.g.,motionvectorsintwodimensions)thereisananalogoustwo-passdistancetransformalgorithmthatcanbeused(e.g.[3]).Messageupdatesusingthetruncatedlinearmodelinequation(4)cannoweasilybecomputedin(k)time.Notethatatruncatedlinearfunctionisthelowerenvelopeofalinearfunctionandtheconstantfunctiondenedbythetruncationvalue.Usingalgebraicpropertiesofminconvolu-tions(see[7])wecancomputeamessageunderthetruncatedlinearmodelintermsofamessageunderthelinearmodelandamessageunderaconstantpenaltymodel.Firstwecomputewhatthemessage,m0,wouldbewiththelinearmodelandthencomputetheelement-wiseminimumofthelinearcostmessageandthevalueusedforthePottscomputation,m(fq)=min(m0(fq);minph(fp)+d):8 Figure2:Theminconvolutionasthelowerenvelopeofkparabolas.3.3QuadraticModelAnothercommonlyusedcostfunctionisthetruncatedquadratic.Inthecaseofaone-dimensionallabelsetthecostgrowsproportionallyto(fp fq)2uptosomelevelandthenbecomesaconstantthereafter.Asintheprevioussubsection,werstconsiderthecasewithouttruncation.Substitutingintothemessageupdateequation(3),thesquaredEuclidean(orquadratic)costupdateisgivenby,mtp!q(fq)=minp c(fp fq)2+h(fp):(6)Analogoustothelinearcase,thiscanbeviewedasthelowerenvelopeofacollectionoffunctions.Eachvalueoffpdenesaconstraintthatisanupward-facingparabolarootedat(fp;h(fp)),andtheoverallminimizationisdenedbythelowerenvelopeoftheseparabolasasshowninFigure2.Ouralgorithmforcomputingthisquadraticminconvolutionhastwosteps.Firstwecomputethelowerenvelopeofthekparabolasjustmentioned.Wethenllinthevaluesofm(fq)bycheckingtheheightofthelowerenvelopeateachgridlocationfq.Notethatthisapproachstartswithsomethingdenedonagrid(thevaluesofh),movestoacombinatorialstructuredenedoverthewholedomain(thelowerenvelopeoftheparabolas)andthenmovesbacktovaluesonthegridbysamplingthelowerenvelope.PseudocodeforthewholeprocedureisshowninAlgorithm1.Themainpartofthealgorithmisthelowerenvelopecomputation.Notethatanytwoparabolas9 deningthelowerenvelopeintersectatexactlyonepoint.Thehorizontalpositionoftheintersec-tionbetweentheparabolacomingfromgridpositionqandtheonefrompis,s=(h(p)+cp2) (h(q)+cq2)2cp 2cq:Ifqpthentheparabolacomingfromqisbelowtheonecomingfromptotheleftoftheintersectionpoints,andaboveittotherightofs.Wecomputethelowerenvelopebysequentiallycomputingthelowerenvelopeoftherstqparabolas,wheretheparabolasareorderedaccordingtotheircorrespondinghorizontalgridloca-tions.Thealgorithmworksbycomputingthecombinatorialstructureofthislowerenvelope.Wekeeptrackofthestructureusingtwoarrays.Thehorizontalgridlocationofthei-thparabolainthelowerenvelopeisstoredinv[i].Therangeinwhichthei-thparabolaofthelowerenvelopeisbelowtheothersisgivenbyz[i]andz[i+1].Thevariablejkeepstrackofthenumberofparabolasinthelowerenvelope.Whenconsideringtheparabolafromq,wenditsintersectionwiththeparabolafromv[j](therightmostparabolainthelowerenvelopecomputedsofar).Therearetwopossiblecases,asillustratedinFigure3.Iftheintersectionisafterz[j],thenthelowerenvelopeismodiedtoindicatethattheparabolafromqisbelowallothersstartingattheintersectionpoint.Iftheintersectionisbeforez[j]thentheparabolafromv[j]shouldnotbepartofthenewlowerenvelope,sowedecreasejtodeletethatparabolaandrepeattheprocedure.Thisalgorithmisasimplerversionofatechniqueforincrementallycomputingtheloweren-velopeofkparabolasin(klogk)time[6].Thatalgorithmoperatesbysortingtheparabolasintoanappropriateordertobeinsertedintothelowerenvelopeinamortizedconstanttime.Inourcasetheproblemissimplerbecausetheparabolasareallofthesameshapeandtheyarealreadysortedintoanappropriateorder.Wenotethatatwo-dimensionalquadraticminconvolutioncanbecomputedbyrstperformingaone-dimensionalminconvolutionalongeachcolumnofthegrid,andthenperformingaone-dimensionalminconvolutionalongeachrowoftheresult(see[7]).Thisargumentextendstoarbitrarydimensions,resultinginthecompositionofone-dimensionalminconvolutionsalongeachdimensionoftheunderlyinggridoflabels.Forstereoandimagerestorationthelabelspaceis10 AlgorithmDT(h)1.j 0(Indexofrightmostparabolainlowerenvelope)2.v[0] 0(Locationsofparabolasinlowerenvelope)3.z[0] 1(Locationsofboundariesbetweenparabolas)4.z[1] +15.forq=1ton 1(Computelowerenvelope)6.s ((h(q)+cq2) (h(v[j])+cv[j]2))=(2cq 2cv[j])7.ifsz[j]8.thenj j 19.goto610.elsej j+111.v[j] q12.z[j] s13.z[j+1] +114.j 015.forq=0ton 1(Fillinvaluesofminconvolution)16.whilez[j+1]q17.j j+118.Dh(q) c(q v[j])2+h(v[j])Algorithm1:TheminconvolutionalgorithmforthesquaredEuclideandistanceinone-dimension.one-dimensional.Inotherearlyvisionproblemssuchasmotionestimationthelabelspaceistwo-dimensional.Asinthelinearcase,messageupdatesusingatruncatedquadraticmodelcanalsobecomputedin(k)time.Againwerstcomputewhatthemessagewouldbewiththequadraticmodelandthencomputetheelement-wiseminimumofthismessagewiththevaluefromthePottscompu-tation.Moreoverwecancomputemessagesunderadiscontinuitycostfunctiondenedbythelowerenvelopeofasmallnumberoflinearandquadraticfunctionsasdescribedin[7].Notealso11 \n \n \r \n Figure3:Thetwopossiblecasesconsideredbythealgorithmwhenaddingtheparabolafromqtothelowerenvelopeconstructedsofar.In(a)s�z[j]whilein(b)sz[j].thatthealgorithmforthequadraticcostfunctioncouldeasilybemodiedtohandleanyconvexdiscontinuitycost.4BPontheGridGraphInthissectionwedescribehowBPcanbeperformedmoreefcientlyforabipartitegraphwhilegettingessentiallythesameresultsasthestandardalgorithm.Thisisanalogoustothered-blacktechniquesusedforGauss-Seidelrelaxations.ThemainissueinusingsuchatechniqueinthecontextofBPisestablishingthatitcomputesthecorrectmessages.Recallthatabipartitegraphisonewherethenodescanbesplitintotwosetssothateveryedgeconnectspairsofnodesindifferentsets.Ifwecolorthegridgraphinacheckerboardpatterneveryedgeconnectsnodesofdifferentcolors,sothegridgraphisbipartite.Themainobservationisthatforabipartitegraphwithnodes[B,whencomputingthemessagesdenedinequation(2)themessagessentfromnodesinonlydependonthemessagessentfromnodesinBandviceversa.Inparticular,ifweknowthemessagessentfromnodesinatiterationt,wecancomputethemessagesfromnodesinBatiterationt+1.Atthispointwecancomputethemessagesfromnodesinatiterationt+2.Thusthemessagesfromnodesinat12 iterationt+2canbecomputedwithoutevercomputingthemessagesfromthosenodesatiterationt+1.ThismotivatesthefollowingmodicationofthestandardBPalgorithmforbipartitegraphs.Inthenewschememessagesareinitializedinthestandardway,butwealternatebetweenupdatingthemessagesfromandthemessagesfromB.Forconcretenessletmtp!qbethemessagesentfromnodeptonodeqatiterationtunderthisnewmessagepassingscheme.WhentisoddweupdatethemessagessentfromnodesinandkeeptheoldvaluesforthemessagessentfromnodesinB.WhentisevenweupdatethemessagessentfromBbutnotthosesentfrom.Soweonlycomputehalfthemessagesineachiteration.Moreoverwecanstorenewmessagesinthesamememoryspacewheretheoldmessageswere.Thisisbecauseineachiterationthemessagesbeingupdateddonotdependoneachother.Usingtheideasfromthelastparagraphitisstraightforwardtoshowbyinductionthatforallt�0,iftisodd(even)thenmtp!q=8�&#x-0.1;è¹â:mtp!qifp2(ifp2B)mt 1p!qotherwise:Thatis,themessagesmsentunderthenewschemearenearlythesameasthemessagesmsentunderthestandardscheme.NotethatwhenBPconverges,thisalternativemessagepassingschemeconvergestothesamexedpoint.Thisisbecauseafterconvergencemt 1p!q=mtp!q.5Multi-GridBPOnedrawbackofusingBPformanyearlyvisionproblemsfollowsfromthefactthatmessagesareupdatedlocallyandinparallel(atleastconceptually,eventhoughtheimplementationisusuallysequential).Thisimpliesthatittakesmanyiterationsforinformationtoowoverlargedistancesinthegridgraph.Inthissectionwedescribeamulti-gridtechniquetocircumventthisproblem.ThebasicideaistoperformBPinacoarse-to-nemanner,sothatlongrangeinteractionsbe-tweenpixelscanbecapturedbyshortpathsincoarsegraphs.WhilehierarchicalBPmethodshavebeenusedinotherworksuchas[14],ourmethoddiffersinthatweusethehierarchyonlytoinitial-izemessagesatsuccessivelynerlevels.Thismakesitpossibletoreducethenumberofmessagepassingiterationsateachlevel,withoutchangingtheoverallproblemstructure.Incontrast,for13 examplein[14]theunderlyinggraphischangedsoastoreplaceedgesbetweenneighboringpix-elsintheimagegridbyedgesbetweenapixelanditsparentinaquad-treestructure.Thishasthenicepropertyofremovingloopsfromthegraph,butitalsosubstantiallychangestheminimizationproblembeingsolved.Inparticular,thequad-treestructurecreatesartifactsduetothespatiallyvaryingneighborhoodstructure.BPworksbylookingforxedpointsofthemessageupdaterule.Formax-productBPthemessagesareusuallyinitializedtozero(inlog-probabilityspace).Ifwecouldinitializethemes-sagesclosetoaxedpointonewouldexpecttogetconvergencemorerapidly.Thisishowthemethoddescribedhereworks;werunBPatonelevelofresolutionandthenusethemessagesatthatlevelinordertogetestimatesforthemessagesatthenextnerlevel.Thusthecoarse-to-necomputationspeedsupconvergenceoftheoriginalBPproblemleavingthegraphstructureandtheenergyfunctionunchanged.Indevelopingthemulti-gridapproachweuseaslightlydifferentnotationthatmakestheimagegrid explicit.Theproblemthatwewanttosolveistoassignalabelfi;j2Ltoeachlocation(i;j)2 whileminimizingtheenergy,E(f)=X(i;j)2 i;j(fi;j)+X(i;j)2 nCV(fi;j fi+1;j)+X(i;j)2 nRV(fi;j fi;j+1);(7)whereCandRarerespectivelythelastcolumnandlastrowoftheimagegrid.Equation(7)isthesameastheoriginalenergyfunctionin(1)exceptthatitisexpressedoverthelocationsintheimagegridratherthanoverasetofsitesandneighbors.Let 0; 1;:::beahierarchyofgridssuchthat 0= andeachnodein `correspondstoablockofpixelsoftheoriginalgrid ,where=2`.Intuitivelythe`-thlevelrepresentslabelingswheretheimagepixelsineachblockareassignedthesamelabel.Akeypropertyofthisconstructionisthatlongrangeinteractionscanbecapturedbyshortpathsinthecoarselevelgrids,asthepathsgothroughblocksinsteadofgoingthroughpixels.Figure4illustratestwolevelsofthestructure.Nowwedeneahierarchyofoptimizationproblemsonthecoarsegrids.Letf`bealabelingforthesitesin `.Theenergyfunctionatlevel`isgivenby,E(f`)=X(i;j)2 ``i;j(f`i;j)+X(i;j)2 `nC`V`(f`i;j f`i+1;j)+X(i;j)2 `nR`V`(f`i;j f`i;j+1);(8)14 Figure4:Illustrationoftwolevelsinthemulti-gridmethod.Eachnodeinlevel`correspondstoa22blockofnodesinlevel` 1.where`andV`arethedataanddiscontinuitycostsatlevel`.Thereareanumberofoptionsforhowtodenethecostsateachlevel.Wetakeanapproachmotivatedbynite-elementmethods,wherethefullsetofimagepixelscorrespondingtoeachblockistakenintoconsideration.Firstconsiderthedatacost`i;j.Intuitivelyassigningalabel toablock(i;j)atlevel`isequivalenttoassigningthatlabeltoeachpixelintheblock,yieldingasumofthedatacostsforthepixelsinthatblock,`ij( )= 1Xu=0 1Xv=0i+u;j+v( ):Thesummationofnegativelogcostscorrespondstotakingaproductofprobabilities,thusthedatacostforanblockcanbeunderstoodintermsoftheprobabilityofobservingthecorrespondingsetofpixelsgivenoneparticularlabelforallofthem.Agivenblockcanpreferseverallabels,becauseacostisdeterminedforeachlabeloftheblock.Forinstance,ifhalfthepixelspreferlabel andhalfpreferlabel ,theneachoftheselabelswillhavelowcostwhereasotherlabelswillhavehighcost.Notethatwhencomputingthedatacostsitisnotnecessarytoalwayssumovertheoriginalgrid .Insteadthecalculationcanbedonemoreefcientlybysummingoverfourdatacostsatthenextnerlevel.Nowconsiderthediscontinuitycostsatlevel`.Thereisnodiscontinuitycostbetweenpixelsinsideablock,aseverycoarselabelingassignsthesamelabelforsuchpixels.Foreachpairofneighboringblockstherearepairsofpixelsalongtheirboundary.Inmeasuringthedifferencebetweenlabelsfortwoneighboringblocksweuseanitedifferenceapproach,wherethedifference15 betweenthelabelsisdividedbytheseparationbetweentheblockcenters,.Thisleadsto,V`( )=V :ThetermmultiplyingVtakesintoaccountthenumberofneighboringpixelsalongtheboundaryoftwoneighboringblocks,whiletheterminthedenominatorinsideVtakesintoaccounttheseparationbetweenblockswhenmeasuringthedifferencebetweenneighboringlabels.Differentformsofdiscontinuitycostsproducedifferentrelationshipsbetweenthediscontinuitycostsacrosstheproblemhierarchy.Forinstance,usingalinearcostfunctionV(x)=cxyieldsahierarchicaldiscontinuitycostthatisindependentofthelevel,V`(x)=cx;asthetermscancelout.OntheotherhandusingaquadraticcostfunctionV(x)=cx2yieldsahierarchicaldiscontinuitycostthatisweakerhigher-upinthehierarchy,V`(x)=cx2=;Asmentionedbefore,inpracticeitisimportanttouserobustdiscontinuitycostssuchasthetruncatedlinearmodelin(4).Wedothisbytruncatingthediscontinuitycostsateachlevel,V`( )=minV ;d:AnotheralternativewouldbetotruncatetheindividualcostfunctionsV,butthiswouldresultinthetruncationfactorchangingbasedonthelevelinthehierarchy,duetothemultiplicationofVby.Inpracticewehavefounditbettertotruncatethecostsbetweenblocksinstead.Asimplecoarse-to-nestrategyusingthehierarchyofproblemsdenedbyequation(8)istocomputetheBPmessagesfortheproblematthecoarsestlevelofthehierarchyandthenusethattoinitializethemessagesatthenextlevel,andsoon,downtotheoriginalgrid.Themessagesateachlevelareafunctionofthesamesetoflabelsbutrepresentinteractionsbetweendifferentsizedblocksofpixels.Givenanalsetofmessagessentbyablockatlevel`,weinitializethemessagessentbythefourblocksinsideitatlevel` 1tothosevalues.Thisisdoneseparatelyforthemessagesinthefourdirections:right,left,upanddowninthegrid.16 Wehavefoundthatwiththiscoarse-to-neapproachitisenoughtorunBPforasmallnumberofiterationsateachlevel(betweenveandten).Notethatthetotalnumberofnodesinthehierarchyisjust4=3thenumberofnodesatthenestlevel.Thusforagivennumberofiterationsthetotalnumberofmessageupdatesinthehierarchicalmethodisjust1=3morethanthenumberofupdatesinthenestlevel.Thishierarchicalmethoddiffersinasubtlebutimportantwayfromothermulti-scaletech-niquescommonlyusedincomputervision,suchastheGaussianpyramid(e.g.,[5]).Typicallysuchtechniqueshavebeenusedforndingdisplacementsbetweenpixelsinpairsofimagesus-ingdifferentialmethods.Thesetechniquesarebasedonreducingtheresolutionoftheimagedata,whereasoursisbasedonreducingonlytheresolutionatwhichthelabelsareestimated.Forinstanceconsidertheproblemofstereo.Reducingtheimageresolutionreducesthenumberofdisparitiesthatcanbedistinguished.Bythefourthlevelofsuchahierarchy,alldisparitiesbe-tween0and16areindistinguishable.Incontrastourmethoddoesnotlowertheimageresolutionbutratheraggregatesdatacostsoverlargerspatialneighborhoods.Thusevenataveryhighlevelofthehierarchy,smalldisparitiesarestillevidentiftheyarepresentoveralargespatialregion.Thisdifferenceiscrucialtosolvingtheproblemathand,becausewewanttobeabletopropagateinformationaboutquantitiessuchasdisparitiesoverlargeareasoftheimageinasmallnumberofmessagepassingiterations.Ingeneral,weneedanumberoflevelsproportionaltolog2oftheimagediameter.IncontrastaGaussianpyramidhasnousefulinformationaboutdisplacementsatlevelshigherthanlog2ofthemaximummagnitudedisplacement(andthisvalueisusuallymuchsmallerthantheimagediameter).6Sum-ProductBeliefPropagationThemax-productBPalgorithmismotivatedbyndingalabelingwithmaximumposteriorproba-bility,orequivalentlywithminimumenergy.Anothercommonformulationisbasedonselectingthemostprobablelabelforeachpixel.Thereisasubtlebutimportantdifferencebetweenselectingthemostprobablelabelingandselectingthemostprobablelabelforeachpixelindividually.Se-lectingthemostprobablelabelforeachpixelminimizesthenumberofpixelswithincorrectlabels,17 buttheoveralllabelingobtainedinthiswaycouldhavesmalljointposteriorprobability.Thesum-productBPalgorithmcanbeusedtoapproximatetheposteriorprobabilityofeachlabelforeachpixel.Aswiththemax-productalgorithm,thesum-productmethodworksbypassingmessagesaroundthegraphdenedbytheneighborhoodsystemN.Inthissectionweletmtp!qbethemessagethatnodepsendstoaneighboringnodeqatiterationtofthesum-productalgorithm,mtp!q(fq)=Xp0@ (fp fq)p(fp)Ys2N(p)nqmt 1s!p(fp)1A(9)whereasaboveN(p)nqdenotestheneighborsofpotherthanq.Thepotentialfunctionsaredenedintermsofthediscontinuitycostsanddatacostsintheenergyfunction(1),p(fp)=e Dp(p); (fp fq)=e (p ):AfterTiterationsabeliefvectoriscomputedforeachnode,q(fq)=q(fq)Yp2N(q)mTp!q(fq):Thevalueq(fq)isanapproximationtotheprobabilitythatthecorrectlabelforpixelqisfq.Aswastrueforthemax-productcase,thestandardimplementationofthismessagepassingalgorithmonthegridgraphrunsin(nk2T)time,wherenisthenumberofpixelsintheimage,kisthenumberofpossiblelabelsforeachpixelandTisthenumberofiterationsAllofthealgorithmictechniquesthatwediscussedaboveformax-productalsoapplytothesum-productalgorithmforlow-levelvisionproblems.ThebipartitegraphtechniqueinSection4andthemulti-gridtechniqueinSection5bothapplydirectly,asneithertechniquedependsontheparticularformofthemessages.Ontheotherhand,thetechniquesforlinear-timemessageupdatesdependontheformofthemessageandthusdonotapplydirectly.Howeverthereisananalogoussetoftechniquesthatwenowdescribe.Followingtheanalysisforthemax-productcase,wecanrewritethemessageupdateruleinequation(9)asmtp!q(fq)=Xp( (fp fq)h(fp));(10)18 whereh(fp)=p(fp)Qmt 1s!p(fp).Inthisformweseethatthemessageupdatecomputationisaconvolution,whichcanbecomputedin(klogk)timeforkdiscretevaluesoffpandfqusingthefastfouriertransform(FFT).Themostcommonlyusedcompatibilityfunctions (fp fq)areGaussiansormixturesofGaussians,andinthesecasesthemessagecomputationcanbeapproximatedin(k)time.ThemethodisalsoveryfastinpracticeandthuspreferabletotheFFTnotonlybecauseofthelogfactorimprovementbutbecauseofthelowconstants.ConvolutionwithaGaussiancanbeapproximatedin(k)timeusingtheboxsummethodin[13].ThetechniqueusesthefactthataGaussiancanbeaccuratelyapproximatedbythesequentialconvolutionofasmallnumberofboxlters.Thediscreteconvolutionofafunctionwithksam-plepointswithaboxltercanbecomputedin(k)time,becauseeachsuccessiveshiftinvolvesonlyasingleadditionandsubtraction,regardlessofthewidthofthebox.ToapproximateGaus-sianconvolutiontheinputfunctionh(fp)issequentiallyconvolvedwithasetofsuchboxlters.InpracticeonlyfourconvolutionsarenecessarytoobtainagoodapproximationtoaGaussian,yieldingan(k)methodthatisalsoveryfastinpractice.Usingthebox-sumtechniquetogetherwiththemulti-gridandbipartitegraphtechniquesresultsinan(nk)algorithmforsum-productbeliefpropagationonagridwithnnodes(orpixels).Formoregeneralpotentialfunctions ,theuseoftheFFTyieldsan(nklogk)method.7ExperimentsInthissectionweshowsomesimpleexperimentalresultstoillustratethetechniquesdescribedinthepaper.Intheseexperimentsweusedthemax-productformulationofbeliefpropagation,ormorepreciselythemin-sumalgorithmwherecostscorrespondtonegativelogprobabilities.Weconsideredboththeproblemsofstereoandimagerestoration.Inbothcaseswecombinedallthreetechniquestogether:thelineartimemessageupdates,thebipartitegraphmessagepassingschedule,andthemulti-gridmethod.Forthemulti-gridmethodweusedsixlevelsinthegridhierarchy.Thetestimagesweregenerallyaround640480pixelsinsize,makingthecoarsestgridhavejustafewblocks.19 Forthestereoproblemthelabelscorrespondtodifferentdisparitiesthatcanbeassignedtopixelsintheleftimage.Weusedatruncatedlinearcostfunctionforthediscontinuityterm,V(fp fq)=min(fp fq;d):Usingthebrightnessconstancyassumptionweexpectthatcorrespondingpixelsintheleftandrightimagesshouldhavesimilarintensities.Weassumethattheimageshavebeenrectiedsothatdisparitiesarehorizontalshifts,anduseatruncatedlinearmodelforthedatacostp(fp)=min(Il(x;y) Ir(x fp;y); );where isatruncationvalueandisascalingfactor.Thetruncationmakesthedatacostrobusttoocclusionandartifactsthatviolatethebrightnessconstancyassumption(suchasspecularities).Thescalingfactorallowsustocontroltherelativeweightingofthedataanddiscontinuitycosts.Moreinvolveddatacostssuchasthosein[1]couldalsobeemployedinthisframework.Figure5showsstereoresultsfortheTsukubaimagepairusingouralgorithmwiththesetrun-catedlinearcostfunctions.Inthiscaseweused10messageupdateiterationsperlevel.Therunningtimewasaboutonesecondona2GHzPentiumIV.Ingeneratingstereoresultsweusedthefollowingsetofparameters:d=1:7, =15and=0:07.TheresultingdiscontinuitycostfunctionisnearlylikeaPottsmodel,withcostzerowhenthelabelsarethesame,1whentheydifferbyone,and1.7otherwise.Theinputimagesweresmoothedslightly,withaGaussianof=0:7priortocomputingthedatacost.Thisexampleillustratesthattheuseofthehierarchi-calmethodseemstoproducelessvariationintheoutputthanisobtainedbynon-hierarchicalBPtechniques(forexamplethebackgroundismoreuniformlyasingledisparity).Figure6showsthevalueoftheenergythatisbeingminimizedasafunctionofthenumberofmessageupdateiterationsforourmulti-gridBPmethodversusthestandardalgorithm.Notehowourmethodcomputesalowenergysolutioninjustafewiterationsperlevel,whilethestandardalgorithmtakesmanymoreiterationstoobtainasimilarresult.Thisprovidessomeempiricalevidencethatthemulti-gridtechniqueisoperatingasintended,allowinginformationtopropagateoverlongdistancesinfewmessageupdateiterations.Figure7givesempiricalresultsofthespeedupsobtainedbyeachofthetechniquesdescribedinthepaper.ThegraphcomparesrunningBPwithallspeeduptechniquesversusrunningBP20 Figure5:StereoresultsfortheTsukubaimagepair.iterationsenergy multiscalestandardFigure6:Energyofstereosolutionasafunctionofthenumberofmessageupdateiterations.withallbutoneofthetechniques.Ineachcasetherunningtimeofthealgorithmiscontrolledbyvaryingthenumberofmessageupdateiterations.Weseethateachspeeduptechniqueprovidesasignicantbenet.Notehowtheminconvolutionmethodprovidesanimportantspeedupevenwhenthenumberoflabelsissmall(16disparitiesfortheTsukubaimages).Table1showsevaluationresultsofourstereoalgorithmontheMiddleburystereobenchmark[9].Theseresultswereobtainedusingtheparametersdescribedabove.Overallourmethodper-21 Figure7:Energyofstereosolutionasafunctionofrunningtime.ThegraphcomparesrunningBPwillallspeeduptechniquesdescribedinthepaperversusBPwithallbutoneofthetechniques.TsukubaSawtoothVenusMapRankErrorRankErrorRankErrorRankError131.84130.9480.94110.36Table1:EvaluationofthestereoalgorithmontheMiddleburyStereobenchmark.Theerrormea-suresthepercentageofpixelswithwrongdisparities.Ourmethodranksin12thplaceintheoverallevaluation.formscomparablytotheoriginalgraphcutsenergyminimizationapproach[4,9]thatsimilarlyusedsimpledataanddiscontinuitycosts,aswellastoresultsin[11]thatcomparedbeliefpropa-gationwithgraphcuts.Howevertheseothermethodsrunconsiderablymoreslowly,takingtensorhundredsoftimeslongerthanouralgorithm.Itisimportanttostressthatthiscomparisonisintendedtodemonstratethatthealgorithmictechniqueswehavepresentedhereproducesimilarqualityresultsmuchfasterthantheseothermethods.Considerablymoresophisticateddataterms,useofocclusioninformation,andothertechniquescouldbeincorporatedinordertoimprovethe22 accuracyofthenalresults.Whilebeliefpropagationandgraphcutsmethodsarenowcommonlyusedtosolvethestereoproblem,thesetechniqueshavenotbeenwidelyadoptedforotherlowlevelvisionproblemssuchasimagerestoration.ThereisalonghistoryofMRF-basedformulationsofimagerestorationproblems(e.g.,see[2,8]),howevercomputingsolutionsfortheseproblemsusingpreviousmeth-odsisquiteslow,particularlywhenthenumberofpossiblelabelsforeachpixelislarge.Hereweconsidertheproblemofrestoringimagesthathave256intensityvalues.WeuseinputimagesthathavebeencorruptedwithadditiveGaussiannoiseaswellasbymaskingoutregions.Inimagerestorationboththelabelsandtheinputdataareintensities.Weusedatruncatedquadraticforthediscontinuitycost,V(fp fq)=min((fp fq)2;d);andaquadraticcostforthedataterm,p(fp)=min((I(p) fp)2);measuringthedifferencebetweenthelabelatapixelandtheobservedintensityatthatpixel.Inprinciplethisformulationoftherestorationproblemshouldalsodoagoodjobofllinginmissingdata,bypropagatinginformationfromotherpartsoftheimage.TodemonstratethiscapabilityweshowanexampleofanimagethatwasdistortedbyaddingGaussiannoiseof=20,andinwhichinadditionarectangularregionwasmaskedout.Amodieddatacostfunctionwasused,whereforamaskedpixelthedatacostiszeroforanylabel.Thatis,p(fp)=0whenpixelpismasked.Thediscontinuitycostfunctionremainedunchanged.Theparametersforthecostfunctionswere:d=200and=0:04.Inthiscaseweused5messagepassingiterationsperlevelandtherunningtimewasapproximately10secondsona2GhzPentiumIV.Figure8showstheresultsofouralgorithm.Notethatmethoddoesagoodjobofllinginmissingdatabasedontheremainingimage.23 CorruptedRestorationFigure8:Restorationresultswithaninputthathasmissingvalues.8SummaryWehavepresentedthreealgorithmictechniquesforspeedingupthebeliefpropagationapproachforsolvinglowlevelvisionproblemsformulatedintermsofMarkovrandomelds.Themainfocusofthepaperisonthemax-productformulationofbeliefpropagation,andthecorrespondingenergyminimizationproblemintermsofcoststhatareproportionaltonegativelogprobabilities.Wealsoshowhowsimilartechniquesapplytothesum-productformulationofbeliefpropagation.Theuseofourtechniquesyieldsresultsofcomparableaccuracytootheralgorithmsbuthundredsoftimesfaster.InthecaseofstereowequantiedthisaccuracyusingtheMiddleburybenchmark.Themethodisquitestraightforwardtoimplementandinmanycasesshouldremovetheneedtochoosebetweenfastlocalmethodsthathaverelativelylowaccuracy,andslowglobalmethodsthathavehighaccuracy.Therstofthethreetechniquesreducesthetimenecessarytocomputeasinglemessageupdatefrom(k2)to(k),wherekisthenumberofpossiblelabelsforeachpixel.Forthemax-productformulationthistechniqueisapplicabletoproblemswherethediscontinuitycostforneighboringlabelsisatruncatedlinearortruncatedquadraticfunctionofthedifferencebetweenthelabels.Themethodisnotanapproximation,itusesanefcientalgorithmtoproduceexactlythesameresultsasthebruteforcequadratictimemethod.Forsum-productasimilartechniqueyieldsan(klogk)24 methodforanydiscontinuitycostfunctionbasedondifferencebetweenlabels.Thesecondofthethreetechniquesusesthefactthatthegridgraphisbipartitetodecreaseboththestoragerequirementsandtherunningtimebyafactoroftwo.Thisisparticularlyimportantbecauseoftherelativelyhighmemoryrequirementsofbeliefpropagationmethods,whichstoremultipledistributionsateachpixel.Thethirdofthetechniquesusesamulti-gridapproachtoreducethenumberofmessagepassingiterationstoasmallconstant,whereasthestandardmethodrequiresanumberofiterationsthatisproportionaltothediameteroftheimagegrid.Forproblemssuchasstereo,wherethelabelsetisrelativelysmall,thetechniquespresentedhereprovidesubstantialspeedup.Forotherproblems,includingimagerestoration,wherethelabelsetcanbequitelarge,thesetechniquescanmakeanMRF-basedapproachtractablewhereitwasnotbefore.Thereareseveralopportunitiesforfurtherdevelopmentofourtechniques.First,ageneralmethodforcomputingtheminconvolutionquickly,analogoustotheFFTforconvolution,wouldbroadentheapplicabilityoffastmessageupdatestoarbitrarydiscontinuitycostfunctionsbasedondifferencebetweenlabels.Second,thelowerenvelopemethodthatwehavepresentedfortheminconvolutioncouldbeextendedtohandleproblemswherethelabelsareembeddedinsomespacebutdonotlieonaregularlyspacedgrid.Moregenerally,itwouldbeinterestingtoconsiderwhetherothersortsofstructuresonthesetoflabelsenablefastmethods.References[1]S.BircheldandC.Tomasi.Apixeldissimilaritymeasurethatisinsensitivetoimagesam-pling.IEEETransactionsonPatternAnalysisandMachineIntelligence,20(4):401406,April1998.[2]A.BlakeandA.Zisserman.VisualReconstruction.MITPress,1987.[3]G.Borgefors.Distancetransformationsindigitalimages.ComputerVision,GraphicsandImageProcessing,34(3):344371,1986.25 [4]Y.Boykov,O.Veksler,andR.Zabih.Fastapproximateenergyminimizationviagraphcuts.IEEETransactionsonPatternAnalysisandMachineIntelligence,23(11):12221239,2001.[5]P.J.BurtandE.H.Adelson.Thelaplacianpyramidasacompactimagecode.IEEETransac-tionsonCommunication,31(4):532540,1983.[6]O.DevillersandMGolin.Incrementalalgorithmsforndingtheconvexhullsofcirclesandthelowerenvelopesofparabolas.Inform.Process.Lett.,56(3):157164,1995.[7]P.F.FelzenszwalbandD.P.Huttenlocher.Distancetransformsofsampledfunctions.Septem-ber2004.CornellComputingandInformationScienceTechnicalReportTR2004-1963.[8]S.GemanandD.Geman.Stochasticrelaxation,gibbsdistributions,andthebayesianrestora-tionofimages.IEEETransactionsonPatternAnalysisandMachineIntelligence,6(6):721741,1984.[9]D.ScharsteinandR.Szeliski.Ataxonomyandevaluationofdensetwo-framestereocorre-spondencealgorithms.InternationalJournalofComputerVision,47(1):742,2002.[10]J.Sun,N.N.Zheng,andH.Y.Shum.Stereomatchingusingbeliefpropagation.IEEETrans-actionsonPatternAnalysisandMachineIntelligence,25(7):787800,2003.[11]M.F.TappenandW.T.Freeman.Comparisonofgraphcutswithbeliefpropagationforstereo,usingidenticalMRFparameters.InIEEEInternationalConferenceonComputerVision,2003.[12]Y.WeissandW.T.Freeman.Ontheoptimalityofsolutionsofthemax-productbeliefpropaga-tionalgorithminarbitrarygraphs.IEEETransactionsonInformationTheory,47(2):723735,2001.[13]W.M.Wells.Efcientsysthesisofgaussianltersbycascadeduniformlters.IEEETrans-actionsonPatternAnalysisandMachineIntelligence,8(2):234239,1986.[14]A.S.Willsky.Multiresolutionmarkovmodelsforsignalandimageprocessing.ProceedingsoftheIEEE,90(8):13961458,2002.26