FastImageDeconvolutionusingHyper-LaplacianPriors - PDF document

FastImageDeconvolutionusingHyper-LaplacianPriors DilipKrishnan,Dept.ofComputerScience,CourantInstitute,NewYorkUniversitydilip@cs.nyu.eduRobFergus,Dept.ofComputerScience,CourantInstitute,NewYorkUniversityfergus@cs.nyu.eduAbstractTheheavy-taileddistributionofgradientsinnaturalsceneshaveproveneffectivepriorsforarangeofproblemssuchasdenoising,deblurringandsuper-resolution.Thesedistributionsarewellmodeledbyahyper-Laplacianp(x)ekjj
,typ-icallywith0:50:8.However,theuseofsparsedistributionsmakestheproblemnon-convexandimpracticallyslowtosolveformulti-megapixelimages.Inthispaperwedescribeadeconvolutionapproachthatisseveralordersofmag-nitudefasterthanexistingtechniquesthatusehyper-Laplacianpriors.Weadoptanalternatingminimizationschemewhereoneofthetwophasesisanon-convexproblemthatisseparableoverpixels.Thisper-pixelsub-problemmaybesolvedwithalookuptable(LUT).Alternatively,fortwospecicvaluesof,1=2and2=3ananalyticsolutioncanbefound,byndingtherootsofacubicandquarticpoly-nomial,respectively.Ourapproach(usingeitherLUTsoranalyticformulae)isabletodeconvolvea1megapixelimageinlessthan3seconds,achievingcom-parablequalitytoexistingmethodssuchasiterativelyreweightedleastsquares(IRLS)thattake20minutes.Furthermore,ourmethodisquitegeneralandcaneasilybeextendedtorelatedimageprocessingproblems,beyondthedeconvolu-tionapplicationdemonstrated.1IntroductionNaturalimagestatisticsareapowerfultoolinimageprocessing,computervisionandcomputationalphotography.Denoising[14],deblurring[3],transparencyseparation[11]andsuper-resolution[20],arealltasksthatareinherentlyill-posed.Priorsbasedonnaturalimagestatisticscanregularizetheseproblemstoyieldhigh-qualityresults.However,digitalcamerasnowhavesensorsthatrecordim-ageswithtensofmegapixels(MP),e.g.thelatestCanonDSLRshaveover20MP.Solvingtheabovetasksforsuchimagesinareasonabletimeframe(i.e.afewminutesorless),posesaseverechallengetoexistingalgorithms.Inthispaperwefocusononeparticularproblem:non-blinddeconvolution,andproposeanalgorithmthatispracticalforverylargeimageswhilestillyieldinghighqualityresults.Numerousdeconvolutionapproachesexist,varyinggreatlyintheirspeedandsophistication.Simplelteringoperationsareveryfastbuttypicallyyieldpoorresults.Mostofthebest-performingap-proachessolvegloballyforthecorrectedimage,encouragingthemarginalstatisticsofasetoflteroutputstomatchthoseofuncorruptedimages,whichactasapriortoregularizetheproblem.Forthesemethods,atrade-offexistsbetweenaccuratelymodelingtheimagestatisticsandbeingabletosolvetheensuingoptimizationproblemefciently.IfthemarginaldistributionsareassumedtobeGaussian,aclosed-formsolutionexistsinthefrequencydomainandFFTscanbeusedtorecovertheimageveryquickly.However,real-worldimagestypicallyhavemarginalsthatarenon-Gaussian,asshowninFig.1,andthustheoutputisoftenofmediocrequality.AcommonapproachistoassumethemarginalshaveaLaplaciandistribution.Thisallowsanumberoffast`1andrelatedTV-normmethods[17,22]tobedeployed,whichgivegoodresultsinareasonabletime.However,studies1 -100 -80 -60 -40 -20 0 20 40 60 80 100 -15 -10 -5 0 Gradientlog2 Probability Empirical Gaussian (a=2) Laplacian (a=1) Hyper-Laplacian (a=0.66) Figure1:Ahyper-Laplacianwithexponent=2=3isabettermodelofimagegradientsthanaLaplacianoraGaussian.Left:Atypicalreal-worldscene.Right:Theempiricaldistributionofgradientsinthescene(blue),alongwithaGaussiant(cyan),aLaplaciant(red)andahyper-Laplacianwith=2=3(green).Notethatthehyper-Laplaciantstheempiricaldistributionclosely,particularlyinthetails.ofreal-worldimageshaveshownthemarginaldistributionshavesignicantlyheaviertailsthanaLaplacian,beingwellmodeledbyahyper-Laplacian[4,10,18].Althoughsuchpriorsgivethebestqualityresults,theyaretypicallyfarslowerthanmethodsthatuseeitherGaussianorLaplacianpri-ors.Thisisadirectconsequenceoftheproblembecomingnon-convexforhyper-Laplacianswith1,meaningthatmanyofthefast`1or`2tricksarenolongerapplicable.Instead,standardoptimizationmethodssuchasconjugategradient(CG)mustbeused.Onevariantthatworkswellinpracticeisiterativelyreweightedleastsquares(IRLS)[19]thatsolvesaseriesofweightedleast-squaresproblemswithCG,eachonean`2approximationtothenon-convexproblematthecurrentpoint.Inbothcases,typicallyhundredsofCGiterationsareneeded,eachinvolvinganexpensiveconvolutionoftheblurkernelwiththecurrentimageestimate.Inthispaperweintroduceanefcientschemefornon-blinddeconvolutionofimagesusingahyper-Laplacianimagepriorfor01.Ouralgorithmusesanalternatingminimizationschemewherethenon-convexpartoftheproblemissolvedinonephase,followedbyaquadraticphasewhichcanbeefcientlysolvedinthefrequencydomainusingFFTs.Wefocusontherstphasewhereateachpixelwearerequiredtosolveanon-convexseparableminimization.Wepresenttwoapproachestosolvingthissub-problem.Therstusesalookuptable(LUT);thesecondisananalyticapproachspecictotwovaluesof.For=1=2theglobalminimacanbedeterminedbyndingtherootsofacubicpolynomialanalytically.Inthe=2=3case,thepolynomialisaquarticwhoserootscanalsobefoundefcientlyinclosed-form.BothIRLSandourapproachsolveaseriesofapproximationstotheoriginalproblem.However,inourmethodeachapproximationissolvedbyalternatingbetweenthetwophasesaboveafewtimes,thusavoidingtheexpensiveCGdescentusedbyIRLS.Thisallowsourschemetooperateseveralordersofmagnitudefaster.Althoughwefocusontheproblemofnon-blinddeconvolution,itwouldbestraightforwardtoadaptouralgorithmtootherrelatedproblems,suchasdenoisingorsuper-resolution.1.1RelatedWorkHyper-Laplacianimagepriorshavebeenusedinarangeofsettings:super-resolution[20],trans-parencyseparation[11]andmotiondeblurring[9].Inworkdirectlyrelevanttoours,Levinetal.[10]andJoshietal.[7]haveappliedthemtonon-blinddeconvolutionproblemsusingIRLStosolveforthedeblurredimage.Othertypesofsparseimagepriorinclude:GaussianScaleMixtures(GSM)[21],whichhavebeenusedforimagedeblurring[3]anddenoising[14]andstudent-Tdistributionsfordenoising[25,16].Withtheexceptionof[14],thesemethodsuseCGandthusareslow.Thealternatingminimizationthatweadoptisacommontechnique,knownashalf-quadraticsplit-ting,originallyproposedbyGemanandcolleagues[5,6].Recently,Wangetal.[22]showedhowitcouldbeusedwithatotal-variation(TV)normtodeconvolveimages.Ourapproachiscloselyre-latedtothiswork:wealsouseahalf-quadraticminimization,buttheper-pixelsub-problemisquitedifferent.WiththeTVnormitcanbesolvedwithastraightforwardshrinkageoperation.Inourwork,asaconsequenceofusingasparseprior,theproblemisnon-convexandsolvingitefcientlyisoneofthemaincontributionsofthispaper.Chartrand[1,2]hasintroducednon-convexcompressivesensing,wheretheusual`1normonthesignaltoberecoveredisreplacedwitha`pquasi-norm,wherep1.Similartoourapproach,asplittingschemeisused,resultinginanon-convexper-pixelsub-problem.Tosolvethis,aHuber2 approximation(see[1])tothequasi-normisused,allowingthederivationofageneralizedshrinkageoperatortosolvethesub-problemefciently.However,thisapproximatestheoriginalsub-problem,unlikeourapproach.2AlgorithmWenowintroducethenon-blinddeconvolutionproblem.xistheoriginaluncorruptedlineargrayscaleimageofNpixels;yisanimagedegradedbyblurand/ornoise,whichweassumetobeproducedbyconvolvingxwithablurkernelkandaddingzeromeanGaussiannoise.Weas-sumethatyandkaregivenandseektoreconstructx.Giventheill-posednatureofthetask,weregularizeusingapenaltyfunctionj:jthatactsontheoutputofasetoflters1;:::;fjappliedtox.Aweightingtermcontrolsthestrengthoftheregularization.Fromaprobabilisticperspec-tive,weseektheMAPestimateofx:p(xjy;k)p(yjx;k)p(x),thersttermbeingaGaussianlikelihoodandsecondbeingthehyper-Laplacianimageprior.Maximizingp(xjy;k)isequivalenttominimizingthecostlogp(xjy;k):minNXi=10@ 2(xky)2iJXj=1j(xj)ij1A(1)whereiisthepixelindex,andisthe2-dimensionalconvolutionoperator.Forsimplicity,weusetworst-orderderivativelters1=[1-1]and2=[1-1]T,althoughadditionalonescaneasilybeadded(e.g.learnedlters[13,16],orhigherorderderivatives).Forbrevity,wedenoteFjix(xj)iforj=1;::;J.Usingthehalf-quadraticpenaltymethod[5,6,22],wenowintroduceauxiliaryvariablesw1iandw2i(togetherdenotedasw)ateachpixelthatallowustomovetheFjixtermsoutsidethej:jexpression,givinganewcostfunction:min;wXi 2(xky)2i 2kF1ixw1ik22kF2ixw2ik22
jw1ijjw2ij(2)whereisaweightthatwewillvaryduringtheoptimization,asdescribedinSection2.3.As!1,thesolutionofEqn.2convergestothatofEqn.1.MinimizingEqn.2foraxedcanbeperformedbyalternatingbetweentwosteps,onewherewesolveforx,givenvaluesofwandvice-versa.Thenovelpartofouralgorithmliesinthewsub-problem,butrstwebrieydescribethexsub-problemanditsstraightforwardsolution.2.1xsub-problemGivenaxedvalueofwfromthepreviousiteration,Eqn.2isquadraticinx.Theoptimalxisthus:F1TF1F2TF2 KTKxF1Tw1F2Tw2 KTy(3)whereKxxk.Assumingcircularboundaryconditions,wecanapply2DFFT'swhichdiago-nalizetheconvolutionmatricesF1;F2;K,enablingustondtheoptimalxdirectly:xF1F(F1)F(w1)+F(F2)F(w2)+(=)F(K)F(y) F(F1)F(F1)+F(F2)F(F2)+(=)F(K)F(K)(4)whereisthecomplexconjugateanddenotescomponent-wisemultiplication.Thedivisionisalsoperformedcomponent-wise.SolvingEqn.4requiresonly3FFT'sateachiterationsincemanyofthetermscanbeprecomputed.Theformofthissub-problemisidenticaltothatof[22].2.2wsub-problemGivenaxedx,ndingtheoptimalwconsistsofsolving2Nindependent1Dproblemsoftheform:w=argminwjwj 2(wv)2(5)wherevFjix.Wenowdescribetwoapproachestondingw.2.2.1LookuptableForaxedvalueof,winEqn.5onlydependsontwovariables,andv,hencecaneasilybetabulatedoff-linetoformalookuptable.WenumericallysolveEqn.5for10;000differentvaluesofvovertherangeencounteredinourproblem(0:6v0:6).Thisisrepeatedfordifferentvalues,namelyintegerpowersofp 2between1and256.AlthoughtheLUTgivesanapproximatesolution,itallowsthewsub-problemtobesolvedveryquicklyforany�0.3 2.2.2AnalyticsolutionForsomespecicvaluesof,itispossibletoderiveexactanalyticalsolutionstothewsub-problem.For=2,thesub-problemisquadraticandthuseasilysolved.If=1,Eqn.5reducestoa1-Dshrinkageoperation[22].Forsomespecialcasesof12,thereexistanalyticsolutions[26].Here,weaddressthemorechallengingcaseof1andwenowdescribeawaytosolveEqn.5fortwospecialcasesof=1=2and=2=3.Fornon-zerow,settingthederivativeofEqn.5w.r.twtozerogives:jwj1sign(w)+(wv)=0(6)For=1=2,thisbecomes,withsuccessivesimplication:jwj1=2sign(w)+2(wv)=0(7)jwj1=42(vw)2(8)w32vw2v2wsign(w)=42=0(9)AtrstsightEqn.9appearstobetwodifferentcubicequationswiththe1=42term,howeverweneedonlyconsideroneoftheseasvisxedandwmustliebetween0andv.Hencewecanreplacesign(w)withsign(v)inEqn.9:w32vw2v2wsign(v)=42=0(10)Forthecase=2=3,usingasimilarderivation,wearriveat:w43vw3+3v2w2v3w8 273=0(11)therebeingnosign(w)termasitconvenientlycancelsinthiscase.Hencew,thesolutionofEqn.5,iseither0orarootofthecubicpolynomialinEqn.10for=1=2,orequivalentlyarootofthequarticpolynomialinEqn.10for=2=3.Althoughitistemptingtotrythesamemanipulationfor=3=4,thisresultsina5thorderpolynomial,whichcanonlybesolvednumerically.Findingtherootsofthecubicandquarticpolynomials:Analyticformulaeexistfortherootsofcubicandquarticpolynomials[23,24]andtheyformthebasisofourapproach,asdetailedinAlgorithms2and3.Inboththecubicandquarticcases,thecomputationalbottleneckisthecuberootoperation.AnalternativewayofndingtherootsofthepolynomialsEqn.10andEqn.11istouseanumericalroot-ndersuchasNewton-Raphson.Inourexperiments,wefoundNewton-RaphsontobeslowerandlessaccuratethaneithertheanalyticmethodortheLUTapproach(see[8]forfutherdetails).Selectingthecorrectroots:Giventherootsofthepolynomial,weneedtodeterminewhichonecorrespondstotheglobalminimaofEqn.5.When=1=2,theresultingcubicequationcanhave:(a)3imaginaryroots;(b)2imaginaryrootsand1realroot,or(c)3realroots.Inthecaseof(a),thejwjtermmeansEqn.5haspositivederivativesaround0andthelackofrealrootsimpliesthederivativeneverbecomesnegative,thusw=0.For(b),weneedtocomparethecostsofthesinglerealrootandw=0,anoperationthatcanbeefcientlyperformedusingEqn.13below.In(c)wehave3realroots.ExaminingEqn.7andEqn.8,weseethatthesquaringoperationintroducesaspuriousrootabovevwhenv�0,andbelowvwhenv0.Thisrootcanbeignored,sincewmustliebetween0andv.ThecostfunctioninEqn.5hasalocalmaximumnear0andalocalminimumbetweenthislocalmaximumandv.Henceofthe2remainingroots,theonefurtherfrom0willhavealowercost.Finally,weneedtocomparethecostofthisrootwiththatofw=0usingEqn.13.Wecanusesimilarargumentsforthe=2=3case.Herewecanpotentiallyhave:(a)4imaginaryroots,(b)2imaginaryand2realroots,or(c)4realroots.In(a),w=0istheonlysolution.For(b),wepickthelargerofthe2realrootsandcomparethecostswithw=0usingEqn.13,similartothecaseof3realrootsforthecubic.Case(c)neveroccurs:thenalquarticpolynomialEqn.11wasderivedwithacubingoperationfromtheanalyticderivative.Thisintroduces2spuriousrootsintothenalsolution,bothofwhichareimaginary,thusonlycases(a)and(b)arepossible.Inboththecubicandquarticcases,weneedanefcientwaytopickbetweenw=0andarealrootthatisbetween0andv.WenowdescribeadirectmechanismfordoingthiswhichdoesnotinvolvetheexpensivecomputationofthecostfunctioninEqn.51.Letrbethenon-zerorealroot.0mustbechosenifithaslowercostinEqn.5.Thisimplies: 1Thisrequiresthecalculationofafractionalpower,whichisslow,particularlyif=2.4 jrj 2(rv)2v2 2sign(r)jrj1 2(r2v)0;r0(12)Sinceweareonlyconsideringrootsofthepolynomial,wecanuseEqn.6toeliminatesign(r)jrj1fromEqn.6andEqn.12,yieldingthecondition:r2v(1) (2);v0(13)sincesign(r)=sign(v).Sowrifrisbetween2v=3andvinthe=1=2caseorbetweenv=2andvinthe=2=3case.Otherwisew=0.Usingthisresult,pickingwcanbeefcientlycoded,e.g.lines12–16ofAlgorithm2.Overall,theanalyticapproachisslowerthantheLUT,butitgivesanexactsolutiontothewsub-problem.2.3SummaryofalgorithmWenowgivetheoverallalgorithmusingaLUTforthewsub-problem.AsoutlinedinAlgorithm1below,weminimizeEqn.2byalternatingthexandwsub-problemsTtimes,beforeincreasingthevalueofandrepeating.Startingwithsomesmallvalue0wescaleitbyafactorIncuntilitexceedssomexedvalueMax.Inpractice,wendthatasingleinneriterationsufces(T=1),althoughmorecansometimesbeneededwhenissmall. Algorithm1Fastimagedeconvolutionusinghyper-Laplacianpriors Require:Blurredimagey,kernelk,regularizationweight,exponent(¿0)Require:regimeparameters:0;Inc;MaxRequire:NumberofinneriterationsT.1:=0,x=y2:PrecomputeconstanttermsinEqn.4.3:whileMaxdo4:iter=05:fori=1toTdo6:Givenx,solveEqn.5forallpixelsusingaLUTtogivew7:Givenw,solveEqn.4togivex8:endfor9:=Inc
10:endwhile11:returnDeconvolvedimagex Aswithanynon-convexoptimizationproblem,itisdifculttoderiveanyguaranteesregardingtheconvergenceofAlgorithm1.However,wecanbesurethattheglobaloptimumofeachsub-problemwillbefound,giventhexedxandwfromthepreviousiteration.Likeothermethodsthatusethisformofalternatingminimization[5,6,22],thereislittletheoreticalguidanceforsettingtheschedule.WendthatthesimpleschemeshowninAlgorithm1workswelltominimizeEqn.2anditsproxyEqn.1.TheexperimentsinSection3showourschemeachievesverysimilarSNRlevelstoIRLS,butatagreatlylowercomputationalcost.3ExperimentsWeevaluatethedeconvolutionperformanceofouralgorithmonimages,comparingthemtonumer-ousothermethods:(i)`2(Gaussian)prioronimagegradients;(ii)Lucy-Richardson[15];(iii)thealgorithmofWangetal.[22]usingatotalvariation(TV)normpriorand(iv)avariantof[22]usingan`1(Laplacian)prior;(v)theIRLSapproachofLevinetal.[10]usingahyper-Laplacianpriorwith=1=2;2=3;4=5.NotethatonlyIRLSandourmethoduseapriorwith1.FortheIRLSscheme,weusedtheimplementationof[10]withdefaultparameters,theonlychangebeingtheremovalofhigherorderderivativelterstoenableadirectcomparisonwithotherapproaches.NotethatIRLSand`2directlyminimizeEqn.1,whileourmethod,andtheTVand`1approachesof[22]minimizethecostinEqn.2,usingT=1;0=1;Inc=2p 2;Max=256.Inourapproach,weuse=1=2and=2=3,andcomparetheperformanceoftheLUTandanalyticmethodsaswell.Allrunswereperformedwithmultithreadingenabled(over4CPUcores).5 Weevaluatethealgorithmsusingasetofblurryimages,createdinthefollowingway.7in-focusgrayscalereal-worldimagesweredownloadedfromtheweb.Theywerethenblurredbyreal-worldcamerashakekernelsfrom[12].1%Gaussiannoisewasadded,followedbyquantizationto255discretevalues.Inanypracticaldeconvolutionsettingtheblurkernelisneverperfectlyknown.Therefore,thekernelpassedtothealgorithmswasaminorperturbationofthetruekernel,tomimickernelestimationerrors.Inexperimentswithnon-perturbedkernels(notshown),theresultsaresimilartothoseinTables3and1butwithslightlyhigherSNRlevels.SeeFig.2foranexampleofakernelfrom[12]anditsperturbedversion.OurevaluationmetricwastheSNRbetweentheoriginalimage^xandthedeconvolvedoutputx,denedas10log10k^(^)k2 k^k2,(^x)beingthemeanof^x.InTable1wecomparethealgorithmson7differentimages,allblurredwiththesame1919kernel.ForeachalgorithmweexhaustivelysearchedoverdifferentregularizationweightstondthevaluethatgavethebestSNRperformance,asreportedinthetable.InTable3weevaluatethealgorithmswiththesame512512imageblurredby8differentkernels(from[12])ofvaryingsize.Again,theoptimalvalueofforeachkernel/algorithmcombinationwaschosenfromarangeofvaluesbasedonSNRperformance.Table2showstherunningtimeofseveralalgorithmsonimagesupto30723072pixels.Figure2showsalarger2727blurbeingdeconvolvedfromtwoexampleimages,comparingtheoutputofdifferentmethods.Thetablesandguresshowourmethodwith=2=3andIRLSwith=4=5yieldinghigherqualityresultsthanothermethods.However,ouralgorithmisaround70to350timesfasterthanIRLSdependingonwhethertheanalyticorLUTmethodisused.Thisspeedupfactorisindependentofimagesize,asshownbyTable2.The`1methodof[22]isthebestoftheothermethods,beingofcomparablespeedtooursbutachievinglowerSNRscores.TheSNRresultsforourmethodarealmostthesamewhetherweuseLUTsoranalyticapproach.Hence,inpractice,theLUTmethodispreferred,sinceitisapproximately5timesfasterthantheanalyticmethodandcanbeusedforanyvalueof. Image IRLS IRLS IRLS Ours Ours # Blurry `2 Lucy TV `1 =1/2 =2/3 =4/5 =1/2 =2/3 1 6.42 14.13 12.54 15.87 16.18 14.61 15.45 16.04 16.05 16.44 2 10.73 17.56 15.15 19.37 19.86 18.43 19.37 20.00 19.78 20.26 3 12.45 19.30 16.68 21.83 22.77 21.53 22.62 22.95 23.26 23.27 4 8.51 16.02 14.27 17.66 18.02 16.34 17.31 17.98 17.70 18.17 5 12.74 16.59 13.28 19.34 20.25 19.12 19.99 20.20 21.28 21.00 6 10.85 15.46 12.00 17.13 17.59 15.59 16.58 17.04 17.79 17.89 7 11.76 17.40 15.22 18.58 18.85 17.08 17.99 18.61 18.58 18.96 Av.SNRgain 6.14 3.67 8.05 8.58 7.03 7.98 8.48 8.71 8.93 Av.Time 79.85 1.55 0.66 0.75 354 354 354 L:1.01 L:1.00 (secs) A:5.27 A:4.08 Table1:ComparisonofSNRsandrunningtimeof9differentmethodsforthedeconvolutionof7576864images,blurredwiththesame1919kernel.L=Lookuptable,A=Analytic.Thebestperformingalgorithmforeachkernelisshowninbold.Ouralgorithmwith=2=3beatsIRLSwith=4=5,aswellasbeingmuchfaster.Onaverage,boththesemethodsoutperform`1,demon-stratingthebenetsofasparseprior. Image `1 IRLS Ours(LUT) Ours(Analytic) size =4/5 =2/3 =2/3 256256 0.24 78.14 0.42 0.7 512512 0.47 256.87 0.55 2.28 10241024 2.34 1281.3 2.78 10.87 20482048 9.34 4935 10.72 44.64 30723072 22.40 - 24.07 100.42 Table2:Run-timesofdifferentmethodsforarangeofimagesizes,usinga1313kernel.OurLUTalgorithmismorethan100timesfasterthantheIRLSmethodof[10].4DiscussionWehavedescribedanimagedeconvolutionschemethatisfast,conceptuallysimpleandyieldshighqualityresults.Ouralgorithmtakesanovelapproachtothenon-convexoptimizationprob-6 Original LSNR=14.89t=0.1 LSNR=18.10t=0.8 BlurredSNR=7.31 Ours =2/3SNR=18.96t=1.2 IRLS =4/5SNR=19.05t=483.9 Original LSNR=11.58t=0.1 LSNR=13.64t=0.8 BlurredSNR=2.64 Ours =2/3SNR=14.15t=1.2 IRLS =4/5SNR=14.28t=482.1 Figure2:Cropsfromtwoimages(#1)beingdeconvolvedby4differentalgorithms,includingoursusinga2727kernel(#7).Inthebottomleftinset,weshowtheoriginalkernelfrom[12](lower)andtheperturbedversionprovidedtothealgorithms(upper),tomaketheproblemmorerealistic.Thisgureisbestviewedonscreen,ratherthaninprint.7 Kernel IRLS IRLS IRLS Ours Ours #/size Blurry `2 Lucy TV `1 =1/2 =2/3 =4/5 =1/2 =2/3 #1:1313 10.69 17.22 14.49 19.21 19.41 17.20 18.22 18.87 19.36 19.66 #2:1515 11.28 16.14 13.81 17.94 18.29 16.17 17.26 18.02 18.14 18.64 #3:1717 8.93 14.94 12.16 16.50 16.86 15.34 16.36 16.99 16.73 17.25 #4:1919 10.13 15.27 12.38 16.83 17.25 15.97 16.98 17.57 17.29 17.67 #5:2121 9.26 16.55 13.60 18.72 18.83 17.23 18.36 18.88 19.11 19.34 #6:2323 7.87 15.40 13.32 17.01 17.42 15.66 16.73 17.40 17.26 17.77 #7:2727 6.76 13.81 11.55 15.42 15.69 14.59 15.68 16.38 15.92 16.29 #8:4141 6.00 12.80 11.19 13.53 13.62 12.68 13.60 14.25 13.73 13.68 Av.SNRgain 6.40 3.95 8.03 8.31 6.74 7.78 8.43 8.33 8.67 Av.Time 57.44 1.22 0.50 0.55 271 271 271 L:0.81 L:0.78 (sec) A:2.15 A:2.23 Table3:ComparisonofSNRsandrunningtimeof9differentmethodsforthedeconvolutionofa512512imageblurredby7differentkernels.L=Lookuptable,A=Analytic.Ouralgorithmbeatsallothermethodsintermsofquality,withtheexceptionofIRLSonthelargestkernelsize.However,ouralgorithmisfarfasterthanIRLS,beingcomparableinspeedtothe`1approach.lemarisingfromtheuseofahyper-Laplacianprior,byusingasplittingapproachthatallowsthenon-convexitytobecomeseparableoverpixels.UsingaLUTtosolvethissub-problemallowsforordersofmagnitudespeedupinthesolutionoverexistingmethods.OurMatlabimplementationisavailableonlineathttp://cs.nyu.edu/˜dilip/wordpress/?page_id=122.Apotentialdrawbacktoourmethod,commontotheTVand`1approachesof[22],isitsuseoffrequencydomainoperationswhichassumecircularboundaryconditions,somethingnotpresentinrealimages.Thesegiverisetoboundaryartifactswhichcanbeovercometosomeextendwithedgetaperingoperations.However,ouralgorithmissuitableforverylargeimageswheretheboundariesareasmallfractionoftheoverallimage.Althoughwefocusondeconvolution,ourschemecanbeadaptedtoarangeofotherproblemswhichrelyonnaturalimagestatistics.Forexample,bysettingk=1thealgorithmcanbeusedtodenoise,orifkisadefocuskernelitcanbeusedforsuper-resolution.Thespeedofferedbyouralgorithmmakesitpracticaltoperformtheseoperationsonthemulti-megapixelimagesfrommoderncameras. Algorithm2:SolveEqn.5for=1=2 Require:Targetvaluev,Weight1:=1062:Computeintermediarytermsm;t1;t2;t33:m=sign(v)24:t1=2v=5:t2=3p mv3+3p p m2+4mv36:t3=v2=t27:Compute3roots,1;r2;r38:1=t1+1(3
1=3)
t2+21=3
t39:2=t1(1p i)(6
1=3)
t2(1+p i)(3
2=3)
t310:3=t1(1+p i)(6
1=3)
t2(1p i)(3
2=3)
t311:Pickglobalminimumfrom(0;r1;r2;r3)12:=[1;r2;r3]13:c1=(abs(imag()))Rootmustbereal14:c2=real()sign(v)&#x-285;&#x.667;(2
abs(v))RootmustobeyboundofEqn.1315:c3=real()sign(v)abs(v)Rootv16:w=max((c1&c2&c3)real()sign(v))sign(v)returnw Algorithm3:SolveEqn.5for=2=3 Require:Targetvaluev,Weight1:=1062:Computeintermediarytermsm;t1;:::;t73:m=8(273)4:t1=
v25:t2=v36:t3=
mv27:t4=t32+p m327+m2v48:t5=3p t49:t6=2(
t1+t5+m=(3
t5))10:t7=p t13+t611:Compute4roots,1;r2;r3;r412:1=3v=4+(t7+p (t1+t6+t2=t7))13:2=3v=4+(t7p (t1+t6+t2=t7))14:3=3v=4+(t7+p (t1+t6t2=t7))15:4=3v=4+(t7p (t1+t6t2=t7))16:Pickglobalminimumfrom(0;r1;r2;r3;r4)17:=[1;r2;r3;r4]18:c1=(abs(imag()))Rootmustbereal19:c2=real()sign(v)&#x-285;&#x.670;(1
abs(v))RootmustobeyboundinEqn.1320:c3=real()sign(v)abs(v)Rootv21:w=max((c1&c2&c3)real()sign(v))sign(v)returnw 8 References[1]R.Chartrand.Fastalgorithmsfornonconvexcompressivesensing:Mrireconstructionfromveryfewdata.InIEEEInternationalSymposiumonBiomedicalImaging(ISBI),2009.[2]R.ChartrandandV.Staneva.Restrictedisometrypropertiesandnonconvexcompressivesens-ing.InverseProblems,24:1–14,2008.[3]R.Fergus,B.Singh,A.Hertzmann,S.T.Roweis,andW.Freeman.Removingcamerashakefromasinglephotograph.ACMTOG(Proc.SIGGRAPH),25:787–794,2006.[4]D.Field.Whatisthegoalofsensorycoding?NeuralComputation,6:559–601,1994.[5]D.GemanandG.Reynolds.Constrainedrestorationandrecoveryofdiscontinuities.PAMI,14(3):367–383,1992.[6]D.GemanandC.Yang.Nonlinearimagerecoverywithhalf-quadraticregularization.PAMI,4:932–946,1995.[7]N.Joshi,L.Zitnick,R.Szeliski,andD.Kriegman.Imagedeblurringanddenoisingusingcolorpriors.InCVPR,2009.[8]D.KrishnanandR.Fergus.Fastimagedeconvolutionusinghyper-laplacianpriors,supple-mentarymaterial.NYUTech.Rep.2009,2009.[9]A.Levin.Blindmotiondeblurringusingimagestatistics.InNIPS,2006.[10]A.Levin,R.Fergus,F.Durand,andW.Freeman.Imageanddepthfromaconventionalcamerawithacodedaperture.ACMTOG(Proc.SIGGRAPH),26(3):70,2007.[11]A.LevinandY.Weiss.Userassistedseparationofreectionsfromasingleimageusingasparsityprior.PAMI,29(9):1647–1654,Sept2007.[12]A.Levin,Y.Weiss,F.Durand,andW.T.Freeman.Understandingandevaluatingblinddecon-volutionalgorithms.InCVPR,2009.[13]S.Osindero,M.Welling,andG.Hinton.Topographicproductmodelsappliedtonaturalscenestatistics.NeuralComputation,1995.[14]J.Portilla,V.Strela,M.J.Wainwright,andE.P.Simoncelli.ImagedenoisingusingascalemixtureofGaussiansinthewaveletdomain.IEEETIP,12(11):1338–1351,November2003.[15]W.Richardson.Bayesian-basediterativemethodofimagerestoration.62:55–59,1972.[16]S.RothandM.J.Black.FieldsofExperts:AFrameworkforLearningImagePriors.InCVPR,volume2,pages860–867,2005.[17]L.Rudin,S.Osher,andE.Fatemi.Nonlineartotalvariationbasednoiseremovalalgorithms.PhysicaD,60:259–268,1992.[18]E.SimoncelliandE.H.Adelson.Noiseremovalviabayesianwaveletcoring.InICIP,pages379–382,1996.[19]C.V.Stewart.Robustparameterestimationincomputervision.SIAMReviews,41(3):513–537,Sept.1999.[20]M.F.Tappen,B.C.Russell,andW.T.Freeman.Exploitingthesparsederivativepriorforsuper-resolutionandimagedemosaicing.InSCTV,2003.[21]M.WainwrightandS.Simoncelli.Scalemixturesofgaussiansandtehstatisticsofnaturalimages.InNIPS,pages855–861,1999.[22]Y.Wang,J.Yang,W.Yin,andY.Zhang.Anewalternatingminimizationalgorithmfortotalvariationimagereconstruction.SIAMJ.ImagingSciences,1(3):248–272,2008.[23]E.W.Weisstein.Cubicformula.http://mathworld.wolfram.com/CubicFormula.html.[24]E.W.Weisstein.Quarticequation.http://mathworld.wolfram.com/QuarticEquation.html.[25]M.Welling,G.Hinton,andS.Osindero.Learningsparsetopographicrepresentationswithproductsofstudent-tdistributions.InNIPS,2002.[26]S.Wright,R.Nowak,andM.Figueredo.Sparsereconstructionbyseparableapproximation.IEEETrans.SignalProcessing,pageToappear,2009.9