当前位置：文档之家 › Isolated prompt photon photoproduction at NLO

Isolated prompt photon photoproduction at NLO

  • 格式：pdf
  • 大小：500.54 KB
  • 文档页数：22

a rXiv:h ep-ph/15121v114May21LAPTH-845/01LPT-Orsay 01-43May 2001Isolated prompt photon photoproduction at NLO M.Fontannaz a ,J.Ph.Guillet b ,G.Heinrich a a Laboratoire de Physique Th´e orique 1LPT,Universit´e de Paris XI,Bˆa timent 210,F-91405Orsay,France b Laboratoire d’Annecy-Le-Vieux de Physique Th´e orique 2LAPTH,Chemin de Bellevue,B.P.110,F-74941Annecy-le-Vieux,France Abstract We present a full next-to-leading order code to calculate the photoproduction of prompt photons.The code is a general purpose program of ”partonic event generator”type with large flexibility.We study the possibility to constrain the photon structure functions and comment on isolation issues.A comparison to ZEUS data is also shown.1IntroductionHigh energy electron-proton scattering at the DESY ep collider HERA is dominated by photopro-duction processes,where the electron is scattered at small angles,emitting a quasireal photon which scatters with the proton.These processes are of special interest since they are sensitive to both the partonic structure of the photon as well as of the proton.In particular,they offer the possibil-ity to constrain the(presently poorly known)gluon distributions in the photon,since in a certain kinematical region the subprocess qg→γq,where the gluon is stemming from a resolved photon, is dominating.Up to now,the experimental errors were too large to discriminate clearly between different sets of gluon distributions in the photon,but a high statistics analysis of the1996-2000 HERA data on prompt photon photoproduction announced by the ZEUS collaboration will shed new light on this issue.The calculation of higher order corrections to the Compton processγq→γq has been initiated some time ago[1]–[6].The most recent calculations for prompt photon photoproduction have been done by Gordon/Vogelsang[6]for isolated prompt photon production,Gordon[7]for photon plus jet production and by the group Krawczyk/Zembrzuski[8]for both the inclusive case and γ+jet.However,all of these calculations contain certain drawbacks.In[6],isolation is implemented by adding a subtraction term evaluated in the collinear approximation to the fully inclusive cross section.The programs of[7]and[8]do not contain the full set of NLO corrections.In[7],those parts where thefinal state photon comes from fragmentation of a hard parton were included only at leading order,arguing that isolation cuts will suppress the fragmentation component in any case to a large extent.Moreover,the box contribution has not been included.In[8],higher order corrections are included only for the case where initial andfinal state photons are both direct.So not only the contributions from fragmentation,but also the case where the initial photon is resolved are included at Born level only.However,the box contribution has been taken into account.The calculation presented in this paper takes into account the full NLO corrections to all four subparts.The corresponding matrix elements already have been calculated and tested in previous works[2,9,10].A major advantage of the present code is also given by the fact that it is constructed as a”partonic event generator”and as such is veryflexible.Various sorts of observables matching a particular experimental analysis can be defined and histogrammed for an event sample generated once and for all.This strategy already has been applied to construct NLO codes forγγproduction (DIPHOX)[11]and one or two jets photoproduction[12].The paper is organized as follows.In section2wefirst describe the theoretical framework and the treatment of the infrared singularities.Then we discuss the implementation of isolation cuts and outline the structure of the code.Section3is devoted to phenomenology.We study the effect of isolation,determine the kinematic region which is most sensitive to the gluon distribution in the photon and illustrate the sensitivity of the cross section to the energy of the incoming photon.We give results for inclusive isolated prompt photon production and compare with a recent analysis of ZEUS data[13],before we come to the conclusions in section4.2Theoretical formalism and description of the methodIn this section the general framework for prompt photon photoproduction will be outlined.We will review the contributing subprocesses,the treatment of infrared singularities and the implementation of isolation cuts.2.1The subprocesses contributing at NLOThe inclusive cross section for ep→γX can symbolically be written as a convolution of the parton densities of the incident particles(resp.fragmentation function for an outgoing parton fragmenting into a photon)with the partonic cross sectionˆσdσep→γX(P p,P e,Pγ)= a,b,c dx e dx p dz F a/e(x e,M)F b/p(x p,M p)dˆσab→cX(x p P p,x e P e,Pγ/z,µ,M,M p,M F)Dγ/c(z,M F)(1) where M,M p are the initial state factorization scales,M F thefinal state factorization scale andµthe renormalization scale.The subprocesses contributing to the partonic reaction ab→cX can be divided into four cat-egories which will be denoted by1.direct direct 2.direct fragmentation 3.resolved direct 4.resolved fragmentation.The cases”direct direct”and”resolved direct”correspond to c=γand Dγ/c(z,M F)=δcγδ(1−z)in(1),that is,the prompt3photon is produced directly in the hard subprocess.The cases with”direct”attributed to the initial state photon correspond to a=γ,with Fγ/e approximated by the Weizs¨a cker-Williams formula for the spectrum of the quasireal photonsf eγ(y)=αemylnQ2max(1−y)y .(2)The”resolved”contributions are characterized by a resolved photon in the initial state where a parton stemming from the photon instead of the photon itself participates in the hard subprocess. In these cases,F a/e(x e,M)is given by a convolution of the Weizs¨a cker-Williams spectrum with the parton distributions in the photon:F a/e(x e,M)= 10dy dxγf eγ(y)F a/γ(xγ,M)δ(xγy−x e)(3)Examples of diagrams contributing at Born level to the four categories above are shown in Figs.1 and2.In the case of the”direct direct”part,only the Compton processγq→γq contributes at leading order,at NLO the O(αs)corrections fromγq→γqg resp.γg→γq¯q and the corresponding virtual corrections contribute.We also included the box contribution(Fig.3)into the”direct direct”part since it is known to be sizeable[4],although it is formally a NNLO contribution.In the”direct fragmentation”part,thefinal state photon comes from the fragmentation of a hard parton participating in the short distance subprocess.From a technical point of view,afinal state quark-photon collinear singularity appears in the calculation of the subprocessγg→γq¯q.At higher orders,final state multiple collinear singularities appear in any subprocess where a high p T parton(quark or gluon)undergoes a cascade of successive collinear splittings ending up with a quark-photon splitting.These singularities are factorized to all orders inαs and absorbed,at some arbitrary fragmentation scale M F,into quark and gluon fragmentation functions to a photon,Dγ/c(z,M2F). When the fragmentation scale M F,chosen of the order of the hard scale of the subprocess,is large compared to any typical hadronic scale∼1GeV,these functions behave roughly asα/αs(M2F). Then a power counting argument tells that these fragmentation contributions are asymptotically of the same order inαs as the Born term.A consistent NLO calculation thus requires the inclusion of the O(αs)corrections to these contributions.Note that the singularity appearing in the processγg→γq¯q when thefinal state photon is emitted by the quark and becomes collinear,is subtracted and absorbed by the fragmentation function at the scale M F,as explained above.Therefore both the”direct direct”and the”direct fragmentation”parts separately depend strongly on M F and the attribution of thefinite terms to either of these parts is scheme dependent.Only in the sum of these parts the M F dependenceflattens as expected.The collinear singularities appearing at NLO if the incident photon splits into a collinear q¯q pair are absorbed into the functions F q/γ(xγ,M)at the factorization scale M.(Analogous for theorderFigure1:Examples of direct direct and direct fragmentation contributions at leadingFigure3:The box contributionproton distribution functions F b/p(x p,M p);we will set M p=M in the following.)Thus,by the same reasoning as above for thefinal state,the”initial direct”and”initial resolved”parts separately show a strong dependence on M which cancels out in the sum.Therefore it has to be stressed that only the sum over all four parts has a physical meaning.Figure5illustrates these cancellation mechanisms.The overall reduction of the scale dependence when going from leading to next-to-leading order can be seen in Fig.4.The scales M F and M have been set equal toµ,andµhas been varied between µ=pγT/2andµ=2pγT.One can see that the NLO cross section is much more stable against scale variations,it varies by less than10%in thisµrange.Figure4:Dependence of the total cross section on scale variations.µ=M=M F is varied between µ=pγT/2andµ=2pγT.Figure5:Cancellation of the leading dependence on the fragmentation scale M F between con-tributions from direct and fragmentationfinal states,and on the factorization scale M between parts with direct and resolved initial state.The results are normalized to the total cross section at M F=M=µ=pγT.2.2Treatment of infrared singularitiesThere are basically two methods to isolate the infrared singularities appearing in the calculation at NLO:The phase space slicing method[15]and the subtraction method[16].The method used here follows the approach of[11,14]which combines these two techniques.We will outline the strategy only shortly,for more details we refer to[11].For a generic reaction1+2→3+4+5,at least two particles of thefinal state,say3and4, have a high p T and are well separated in phase space,while the last one,say5,can be soft,and/or collinear to either of the four others.In order to extract these singularities,the phase space is cut into two regions:–part I where the norm p T5of the transverse momentum of particle5is required to be less than some arbitrary value p T m taken to be small compared to the other transverse momenta.This cylinder contains the infrared and the initial state collinear singularities.It also containsa small fraction of thefinal state collinear singularities.–parts II a(b)where the transverse momentum vector of the particle5is required to have a norm larger than p T m,and to belong to a cone C3(C4)about the direction of particle3(4), defined by(η5−ηi)2+(φ5−φi)2≤R2th(i=3,4),with R th some small arbitrary number.C3(C4)contains thefinal state collinear singularities appearing when5is collinear to3(4).–part II c where p T5is required to have a norm larger than p T m,and to belong to neither of the two cones C3,C4.This slice yields no divergence,and can thus be treated directly in4 dimensions.The contributions from regions I and IIa,b are calculated analytically in d=4−2ǫdimensions and then combined with the corresponding virtual corrections such that the infrared singularities cancel, except for the initial(resp.final)state collinear singularities,which are factorized and absorbed into the parton distribution(resp.fragmentation)functions.After the cancellation,thefinite remainders of the soft and collinear contributions in parts I and II a,b,c separately depend on large logarithms ln p T m,ln2p T m and ln R th.When combining the different parts,the following cancellations of the p T m and R th dependences occur:In part I,thefinite terms are approximated by collecting all the terms depending logarithmically on p T m and neglecting the terms proportional to powers of p T m.On the contrary,the R th dependence in the conical parts II a and II b,is kept exactly.This means that an exact cancellation of the dependence on the unphysical parameter R th between part II c and parts II a,b occurs,whereas the cancellation of the unphysical parameter p T m between parts II c,II a,b and part I is only approximate. The parameter p T m must be chosen small enough with respect to pγT in order that the neglected terms can be safely dropped out.On the other hand,it cannot be chosen too small since otherwise numerical instabilities occur.We have investigated the stability of the cross section by varying p T m and R th between0.005and0.1(see Figure6)and accordingly chosen the optimal values p T m=0.05GeV,R th=0.05.Figure6:Dependence of the total cross section on variations of the slicing parameter p T m.Rσdenotes the total(nonisolated)cross section normalized to the total cross section evaluated with p T m=0.005,Rσ(p T m)=σtot(p T m)/σtot(p T m=0.005).One can see that there is a plateau where the cross section is fairly insensitive to variations of p T m.The same study has been made for the dependence on R th,but there the cross section is completely stable within the numerical errors since the R th dependence has been kept exactly in all parts of the matrix element.2.3Implementation of isolation cutsIn order to single out the prompt photon events from the huge background of secondary photons produced by the decays ofπ0,η,ωmesons,isolation cuts have to be imposed on the photon signals in the experiment.A commonly used isolation criterion is the following4:A photon is isolated if, inside a cone centered around the photon direction in the rapidity and azimuthal angle plane,theamount of hadronic transverse energy E hadT deposited is smaller than some value E T maxfixed bythe experiment:(η−ηγ)2+(φ−φγ)2≤R2expE hadT≤E T max (4) Following the conventions of the ZEUS collaboration,we used E T max=ǫpγT withǫ=0.1and R exp=1.Isolation not only reduces the background from secondary photons,but also substantially reduces the fragmentation components,as will be illustrated in section3.1.Furthermore,it is important to note that the isolation parameters must be carefullyfixed in order to allow a comparison between data and perturbative QCD calculations.Indeed a part of the hadronic energy measured in the cone may come from the underlying event;therefore even the direct contribution can be cut by the isolation condition if the latter is too stringent.Let us estimate the importance of this effect and assume that the underlying event one-particle inclusive distribution is given bydn(1)2π4p T ,(5)n(1)being normalized to¯n particles per unit of rapidity.The probability that the isolation condition is fulfilled by a particle from an underlying event isn(1)isol= cone dφdη ∞E T max p T dp T dn(1)2 1+2E T max p T (6) With the ZEUS isolation parameters,E T max=0.5GeV for a photon of pγT=ing¯n=3 and p T ≈0.35GeV extracted from[18],one obtainsn(1)isol≈0.33This estimation is very rough and underestimates the true effect because there is also a non-negligible probability to fulfill the isolation condition with two underlying particles falling into the cone.Only a detailed Monte Carlo description of the underlying events can allow a reliable estimate of this non-perturbative effect.Here we just note that the cut put by ZEUS(E T max≈0.5GeV)is likely to be too low to eliminate any underlying event contamination and therefore makes a comparison between the partonic level QCD predictions and the(hadronic level)data difficult.2.4Features of the codeThe code consists of four subparts corresponding to each of the four categories of subprocesses. For each category,the functions corresponding to the parts I,II a,b,c described in section2.2are integrated separately with the numerical integration package BASES[23].Based on the grid produced by this integration,partonic events are generated with SPRING[23]and stored into an NTUPLE or histogrammed directly.It has to be emphasized that we generatefinal state partonic configurations.Hence this type of program does not provide an exclusive portrait offinal states as given by hadronic event generators like PYTHIA or HERWIG.On the other hand,the latter are only of some improved leading logarithmic accuracy.The information stored in the NTUPLE are the4-momenta of the outgoing particles,their types(i.e.quark,gluon or photon),the energy of the incident photon and, in the fragmentation cases,the longitudinal fragmentation variable associated with the photon from fragmentation.Furthermore a label is stored that allows to identify the origin of the event,e.g.if it comes from a2→2or a2→3process.Based on the information contained in these NTUPLES, suitable observables can be defined and different jet algorithms can be studied.The isolation cuts are included already at the integration level,but the user of the program can turn isolation on or offand vary the input parameters for the isolation cut at will.3Numerical results and comparison to ZEUS dataIn this section we present some numerical results for isolated prompt photon production.We restrict ourselves to the inclusive case,photon+jet production will be discussed in detail in a forthcoming publication.For the parton distributions in the proton we take the MRST2[19]parametrization.Our de-fault choice for the photon distribution functions is AFG[20],for comparisons we also used the GRV[21]distributions transformed to the(7)2E ewhere the sum is over all calorimeter cells,E is the energy deposited in the cell and p z=E cosθ.In order to obtain the”true”photon energy y,corrections for detector effects and energy calibration have to be applied to y JB.These corrections are assumed to be uniform over the whole y range and enter into the experimental systematic error.However,as the background varies with the photon energy y,these corrections may not be uniform.It has to be emphasized that the cross section is very sensitive to a variation of the energy range of the photon.(See Figure12and discussion below.)3.1Numerical results for inclusive prompt photon productionIf not stated otherwise,all plots showing the photon rapidity(ηγ)dependence are integrated over 5GeV<pγT<10GeV and0.2<y=Eγ/E e<0.9.Figure7shows a comparison of the NLO to the leading order result for the isolated cross section dσ/dηγThe importance of the box contribution is clearly visible.The higher order corrections enhance the isolated cross section by about40%.Fig.8shows the rapidity distribution of the full cross section before and after isolation.As already mentioned in section2.3,we used the isolation cuts E T max=ǫpγT withǫ=0.1and R exp =1to match those of the ZEUS collaboration.Fig.8also shows the effect of isolation on thefragmentation part5separately.Isolation reduces the fragmentation component to about6%of the total isolated cross section.In Fig.9the relative magnitude of all four components contributing to dσep→γX/dηγbefore and after isolation is shown.Note that isolation increases the contributions with a direct photon in the final state slightly since there the cut mainly acts on a negative term,which is the one where parton 5is collinear to the photon.It should be emphasized that Figure9has to be read with care since the individual parts have no physical meaning and are very sensitive to scale changes.Nevertheless the dominance of the resolved direct part remains if we choose e.g.µ=M=M F=pγT/2or2pγT.Figure10shows the relative magnitude of contributions from resolved and direct photons in the initial state to the isolated cross section.From the pγT distribution one can conclude that theresolved part dominates the cross section for small values of pγT such that it would be useful to lookat the photon rapidity distribution at pγT=5GeV in order to discriminate between different parton distribution functions in the photon.Since the gluon distribution in the photon is of particular interest,the sensitivity to the gluon in the photon is investigated in Fig.11.One can see that the gluon distribution in the photon starts to become sizeable only for photon rapiditiesηγ>1and dominates over the quark distribution for aboutηγ>2.5.Therefore the region of large photon rapidities and small photon p T is the one where the sensitivity to the gluon in the photon is largest.In order to test further the sensitivity to the gluon,we increased the gluon distribution in the photon uniformly by20%.As can be anticipated from Fig.11,the effect becomes sizeable only forηγ>2and leads to an increase of the cross section by about10%forηγ>2.5.We conclude that in the regionηγ<1,there is basically no sensitivity to the gluon in the photon.However,investigating the photon+jet cross section instead of the inclusive case offers larger possibilities to constrain the gluon in the photon since there one can vary the photon and the jet rapidities in order to single out a kinematic region where the sensitivity is large[24].Figure12shows the effect of a ten percent uncertainty in the”true”bounds of the photon energy y.One can see that a change of the lower bound on y has a large effect,in particular at large photon rapidities.This comes from the fact that the Weizs¨a cker-Williams distribution is large and steeply falling at small y.Increasing the lower bound on y therefore removes a large fraction of the direct events with lower energy initial photons.(y=x e for the direct events and largeηγcorrespond to small x e.)At large photon rapidities the spread due to the use of different parton distribution functions for the photon is smaller than the one caused by a10%variation of the lower bound on y.On the other hand,the region of large photon rapidities is of special interest since there the gluon in the photon is dominating.Therefore a small experimental error in the reconstruction of the ”true”photon energy is crucial in order to be able to discriminate between different sets of parton distribution functions in the photon.It has been tested that the effect of using different proton distribution functions–for example the CTEQ4M or the MRST1set of proton distribution functions–is of the order of3%at most. In all photon rapidity bins the spread is smaller than the one caused by different sets of photon distribution functions(which is about10%at small photon rapidities,see e.g.Fig.14).Thus a discrimination between different sets of photon distribution functions should be possible with the forthcoming full1996-2000data set analysis,where the errors on the data are expected to be small enough.Figure7:Comparison of NLO to LO result for the photon rapidity distribution.Figure8:Effect of isolation on the photon rapidity distribution dσep→γX/dηγfor the full cross section and for the fragmentation components separately.Isolation withǫ=0.1,R exp=1.Figure9:Relative magnitude of all four components contributing to dσep→γX/dηγfor the scale.choiceµ=M=M F=pγTFigure10:Comparison of contributions from resolved and direct photons in the initial state for the photon rapidity and transverse momentum distribution,with isolation.Figure11:Ratio of the contribution from quark resp.gluon distributions in the photon to the fullresolved part.Figure12:Photon rapidity distribution dσep→γX/dηγfor isolated prompt photons integrated over5GeV<pγT <10GeV and different lower bounds on y.Solid line:0.2<y<0.9with AFG photonstructure functions,dotted line:bounds on y changed by about10%,dashed line:0.2<y<0.9 with GRV photon structure functions3.2Comparison with ZEUS dataIn this section we compare our results to the ZEUS1996-97data on inclusive prompt photon photoproduction[13].Figures13and14show the photon p T and rapidity distributions with AFG resp.GRV sets of structure functions for the photon.For the p T distribution the agreement between data and theory is quite good.In the rapidity distribution(Fig.14)the datafluctuate a lot,such that the agreement is still satisfactory.However,it seems that theory underpredicts the data in the backward region,whereas the theoretical prediction tends to be higher at large photon rapidities. The curves of Gordon[7]and Krawczyk/Zembrzuski[8]given in[13]also show this trend.At high ηγthe reason for the difference could be that the isolation cut in the experiment removes more events than in the theoretical(parton level)simulation,as discussed in section2.3.Figure15shows that the discrepancy between theory and data at lowηγcomes mainly from the domain of small photon energies,whereas the discrepancy at largeηγis only present in the range of large photon energies.Note that at largeηγand large y the resolved part dominates and the underlying event could have a large multiplicity.Therefore the isolation criterion could also cut on the non-fragmentation contributions as discussed in section2.3.Figure13:Comparison to ZEUS data of photon p T distribution dσep→γX/dpγT for isolated prompt photons.Results for two different sets of parton distributions in the photon are shown.Figure14:Comparison to ZEUS data of photon rapidity distribution dσep→γX/dηγfor isolated prompt photons.Figure15:Photon rapidity distribution dσep→γX/dηγintegrated over5GeV<pγT<10GeV and different subdivisions of photon energies:(a)0.2<y<0.32,(b)0.32<y<0.5,(c)0.5<y<0.9.4ConclusionsWe have presented a program for prompt photon photoproduction which includes the full next-to-leading order corrections to all contributing subparts.It is a general purpose code of partonic event generator type and as such veryflexible.We used it to study the possibility to constrain the quark and gluon distributions in the photon. It turned out that the sensitivity to the gluon distribution in the photon is negligible in the rapidity range−0.7<ηγ<0.9studied by ZEUS.A discrimination between the AFG/GRV sets of parton distributions in the photon is not possible with the present experimental errors on the ZEUS1996/97data.However,a forthcoming analysis of all1996-2000data announced by the ZEUS collaboration will drastically improve this situation.We have shown that the cross section is very sensitive to small variations of the photon en-ergy range.Therefore a good control of the experimental error on the photon energy fraction y (reconstructed experimentally from the Jacquet-Blondel variable y JB)will be crucial for future comparisons.Despite the largefluctuations of the data,one can say that there is a trend that theory overpre-dicts the data in the forward region.The reason might be that the isolation cut imposed at partonic level in the perturbative QCD calculation does not have the same effect as the experimental one.If the experimental cut is too stringent,a large fraction of the hadronic energy in the isolation cone may come from underlying events,such that experimentally a larger number of events is rejected. We gave a rough estimate of the underlying events to be expected in the isolation cone.The possibilities offered by the study of photon+jet photoproduction will be investigated in a forthcoming publication[24].AcknowledgementsWe would like to thank P.Bussey from the ZEUS collaboration for helpful discussions.G.H. would like to thank the LAPTH for its continuous hospitality.This work was supported by the EU Fourth Training Programme”Training and Mobility of Researchers”,network”Quantum Chromo-dynamics and the Deep Structure of Elementary Particles”,contract FMRX–CT98–0194(DG12-MIHT).References[1]D.W.Duke and J.F.Owens,Phys.Rev.D26,1600(1982).[2]P.Aurenche,A.Douiri,R.Baier,M.Fontannaz and D.Schiff,Z.Phys.C24,309(1984).[3]A.C.Bawa,M.Krawczyk and W.J.Stirling,Z.Phys.C50,293(1991).[4]P.Aurenche,P.Chiappetta,M.Fontannaz,J.Ph.Guillet and E.Pilon,Z.Phys.C56,589(1992).[5]L.E.Gordon and J.K.Storrow,Z.Phys.C63,581(1994).[6]L.E.Gordon and W.Vogelsang,Phys.Rev.D52,58(1995).[7]L.E.Gordon,Phys.Rev.D57,235(1998).[8]M.Krawczyk and A.Zembrzuski,hep-ph/9810253.。

合集下载

A New Antialiasing Approach for Image Compositing

A New Antialiasing Approach for Image Compositing


fJean-Pierre.Braquelaire|Anne.Vialardg@labri.u-bordeaux.fr
Abstract.
In this paper we describe a new antialiasing technique for discrete regions. We propose to use a new model of contours called Euclidean paths providing a construction of an approximation of the underlying real boundary of a discrete region. Euclidean paths are used to compute the blending coe cients of the frontier pixels. This method allows us to superimpose a discrete region on a background while avoiding the well-known staircase e ect.
P p
Integer point Discrete contour Euclidean point associated with P Cell 2: Discrete point and an associated Euclidean point. 2
Let C be a discrete contour. It can be de ned as a sequence of discrete points C = (P 0; P 1 ; : : :; P n ). The sequence of real points C = (p0 ; p1 ; : : :; pn ) where pi is a Euclidean point associated with P i is called Euclidean path associated with C . It is obvious that the sampling of any Euclidean path C associated with a discrete contour C is exactly C . Let us now brie y describe a method of construction that provides a smooth Euclidean path from a discrete contour. As any discrete contour can locally be considered as a discrete line segment, one can approximate the real tangent at each point of the contour path. To obtain the associated Euclidean point, we project the contour point on the local tangent. Repeating this operation along the discrete contour provides an associated Euclidean path which satis es the smoothness requirement. The estimation of local tangents is based on the following de nition of an 8-connected discrete line. De nition 1 8] An 8-connected discrete line D is the set of points (x; y ) of 6 Z2 which satisfy the double inequation Ax ? By < + sup(jAj; jBj) where A; B; 2 6 Z . A The rational B is the slope of D and the integer describes its position in the discrete plane. Intuitively the discrete line is the set of integer points enclosed between the real lines of equation ( ? y = Ax? ?supBjAj;jBj)+1 and y = AxB . Let us consider an 8-connected discrete line D and a given point P of D. Let (A; B; ) be the characteristics of D calculated in the local coordinates system of origin P . Without any loss of generality, we suppose here 0 A B (D is located in the rst octant). The other cases can be inferred by using symmetries. The equation of D is now Ax ? By < + B. Let p be the real point of coordinates (x; y ) = (0; ) with = 1?2 B?B . It can be shown 2 1 that veri es ? 1 < < 2 . In other words, p is a Euclidean point associated with P . One can 2 A consider p as a projection of P on the real line D of equation y = B x + 1?2 B?B . Moreover the 2 line D is at equal distance from both real lines bounding D. This construction of a Euclidean point associated with a point of a discrete line is illustrated on Figure 3.
图像去噪 英文文献及翻译

图像去噪 英文文献及翻译

New Method for Image Denoising while Keeping Edge InformationEdge information is the most important high- frequency information of an image, so we should try to maintain more edge information while denoising. In order to preserve image details as well as canceling image noise, we present a new image denoising method：image denoising based on edge detection。

Before denoising, image's edges are first detected, and then the noised image is divided into two parts: edge part and smooth part. We can therefore set high denoising threshold to smooth part of the image and low denoising threshold to edge part。

The theoretical analyses and experimental results presented in this paper show that，compared to commonly—used wavelet threshold denoising methods, the proposed algorithm could not only keep edge information of an image，but also could improve signal-to—noise ratio of the denoised image.In the wavelet domain，the denoising algorithm based on the threshold filter is widely used, because it’s comparatively effic ient and easy to realize. We can select a threshold according to the characteristic of the image, modifying all of the discrete detail coefficients so as to reduce the noise。

1-9 ∧(1520) and ∑(1385) Photoproductions Based on the New CLAS Data

1-9 ∧(1520) and ∑(1385) Photoproductions Based on the New CLAS Data

·12·IMP&HIRFL Annual Report20141-9Λ(1520)andΣ(1385)Photoproductions Basedon the New CLAS Data∗He JunThe study of nucleon resonance is an important topic in hadron physics.The most important way to study the nucleon resonances is the collision of the photon and nucleon.Up to now,the knowledge about the nucleon resonances around2GeV is still scarce.Recently,the CLAS Collaboration at Jefferson National Accelerator Facility released their exclusive photoproduction cross sections for theΛ(1520),Σ0(1385)andΛ(1405)for energies from near threshold up to a center of mass energy W of2.85GeV with large range of the K production angle[1].Since the threshold for the photoproduction ofΛ(1520)is about2.01GeV,the new experimental data with high precision released by the CLAS Collaboration provide an opportunity to study the nucleon resonances above2GeV.The new experimental data of theΛ(1520)andΣ(1385)photoproductions are analyzed in a Regge-plus-resonance approach combined with the input from the constituent quark model[3].As shown in Figs.1and2,the differential cross sections for both photoproductions are well reproduced in our model.The result suggests that the contact terms are dominant in the photoproductions ofΛ(1520)andΣ(1385).The Regge trajectory is found to be essential for reproducing the experimental data at forward angles,especially at high energies.To describe the behavior of the differential cross section at backward angles,theΛintermediate u channel contribution should be included.Fig.1(color online)The differential cross section dσ/dcosθwith variation of cosθ.The full(black),dashed (red),dotted(blue),dash-dot-dotted(brown),dash-dashed(magenta)and dash-dotted(darkgreen)linesare for full model,contact term,K exchange t channel,u channel,K∗exchange t channel and N(2120).Thefull circle(red)and open circle(blue)are for CLAS13data[1]and LEPS10data[2].∗Foundation item:Major State Basic Research Development Program of China(2014CB845405),National Natural Science Founda-tion of China(11275235,11035006)and Knowledge Innovation Project of Chinese Academy of Sciences(KJCX2-EW-N01).2014IMP&HIRFL Annual Report·13·Fig.2(color online)The differential cross section dσ/dcosθfor theΣ(1385)photoproduction from proton as a function of cosθ.The full(black),dashed(red),dash-dotted(brown),dotted(blue)and dash-dot-dotted(green)lines are for the full model,the contact term,the u channel,the t channel and∆(2000),respectively.The data are from[1].The contributions of nucleon resonances are determined by the radiative and strong decay amplitudes predicted from the constituent quark model[4,5].The results show that the contributions from nucleon resonances are small compared with the contact term,u and t channel contributions,but essential to reproduce the experimental data.In theΣ(1835)photoproduction,the contribution from∆(2000)is much more important than the contributions from other nucleon resonances,and the N(2120)is dominant in theΛ(1520)photoproduction.The data indicate that the∆(2000)and N(2120)are the third∆state with spin-parity5/2+and the third nucleon state with spin-parity 3/2−in constituent quark model,respectively.Other nucleon resonances are found to give negligible contributions in the channel considered in this work.The results are helpful to understand the internal structure of the nucleon resonance around2GeV and the“missing resonance”problem.References[1]K.Moriya,et al.,Phys.Rev.C,88,045201(2013)[2]H.Kohri,H.Fujimura,K.Imai,et al.,Phys.Rev.Lett,104(2010)172001.[3]J.He,X.R.Chen,Phys.Rev.C,86(2012)035204.[4]S.Capstick,W.Roberts,Phys.Rev.D,58(1998)074011.[5]S.Capstick,Phys.Rev.D,46(1992)2864.。

文献翻译-使用不变特征的全景图像自动拼接

文献翻译-使用不变特征的全景图像自动拼接

本科毕业论文英文翻译专业名称学生姓名指导教师毕业时间附录1外文原文Automatic Panoramic Image Stitching using Invariant FeaturesMatthew Brown and David G. Lowe{mbrown|lowe}@cs.ubc.caDepartment of Computer Science,University of British Columbia,Vancouver, Canada.AbstractThis paper concerns the problem of fully automated panoramic image stitching. Though the 1D problem (single axis of rotation) is well studied, 2D or multi-row stitching is more difficult. Previous approaches have used human input or restrictions on the image sequence in order to establish matching images. In this work, we formulate stitching as a multi-image matching problem, and use invariant local features to find matches between all of the images. Because of this our method is insensitive to the ordering, orientation, scale and illumination of the input images. It is also insensitive to noise images that are not part of a panorama, and can recognise multiple panoramas in an unordered image dataset. In addition to providing more detail, this paper extends our previous work in the area by introducing gain compensation and automatic straightening steps.1. IntroductionPanoramic image stitching has an extensive research literature and several commercial applications. The basic geometry of the problem is well understood, and consists of estimating a 3×3 camera matrix or homography for each image. This estimation process needs an initialisation, which is typically provided by user input to approximately align the images, or a fixed image ordering. For example, the PhotoStitch software bundled with Canon digital cameras requires a horizontal or vertical sweep, or a square matrix of images. REALVIZ Stitcher version 4 has a user interface to roughly position the images with a mouse, before automatic registration proceeds. Our work is novel in that we require no such initialisation tobe provided.In the research literature methods for automatic image alignment and stitching fall broadly into two categories—direct and feature based. Direct methods have the advantage that they use all of the available image data and hence can provide very accurate registration, but they require a close initialisation. Feature based registration does not require initialisation, but traditional feature matching methods (e.g., correlation of image patches around Harris corners) lack the invariance properties needed to enable reliable matching of arbitrary panoramic image sequences.In this paper we describe an invariant feature based approach to fully automatic panoramic image stitching. This has several advantages over previous approaches. Firstly, our use of invariant features enables reliable matching of panoramic image sequences despite rotation, zoom and illumination change in the input images. Secondly, by viewing image stitching as a multi-image matching problem, we can automatically discover the matching relationships between the images, and recognise panoramas in unordered datasets. Thirdly, we generate high-quality results using multi-band blending to render seamless output panoramas. This paper extends our earlier work in the area by introducing gain compensation and automatic straightening steps. We also describe an efficient bundle adjustment implementation and show how to perform multi-band blending for multiple overlapping images with any number of bands.The remainder of the paper is structured as follows. Section 2 develops the geometry of the problem and motivates our choice of invariant features. Section 3 describes our image matching methodology (RANSAC) and a probabilistic model for image match verification. In section 4 we describe our image alignment algorithm (bundle adjustment) which jointly optimises the parameters of each camera. Sections 5-7 describe the rendering pipeline including automatic straightening, gain compensation and multi-band blending. In section 9 we present conclusions and ideas for future work.2. Feature MatchingThe first step in the panoramic recognition algorithm is to extract and match SIFT features between all of the images. SIFT features are located at scale-space maxima/minimaof a difference of Gaussian function. At each feature location, a characteristic scale and orientation is established. This gives a similarity-invariant frame in which to make measurements. Although simply sampling intensity values in this frame would be similarity invariant, the invariant descriptor is actually computed by accumulating local gradients in orientation histograms. This allows edges to shift slightly without altering the descriptor vector, giving some robustness to affine change. This spatial accumulation is also important for shift invariance, since the interest point locations are typically only accurate in the 0-3 pixel range. Illumination invariance is achieved by using gradients (which eliminates bias) and normalising the descriptor vector (which eliminates gain).Since SIFT features are invariant under rotation and scale changes, our system can handle images with varying orientation and zoom (see figure 8). Note that this would not be possible using traditional feature matching techniques such as correlation of image patches around Harris corners. Ordinary (translational) correlation is not invariant under rotation, and Harris corners are not invariant to changes in scale.Assuming that the camera rotates about its optical centre, the group of transformations the images may undergo is a special group of homographies. We parameterise each camera by a rotation vector []321θθθθ，，= and focal length f . This gives pairwise homographies j ij i u H u ~~= where1-=j Tj i i ij K R R K H (1) and j i u u ~~, are the homogeneous image positions( []1,~i i i u s u =，where i u is the 2-dimensional image position). The 4 parameter camera model is defined by⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=1000000iii f f K (2)and (using the exponential representation for rotations)[],⨯=i e R i θ[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---=⨯000121323i i i i i i i θθθθθθθ (3)Ideally one would use image features that are invariant under this group of transformations. However, for small changes in image positionj u j i i i u u u u u i ∆∂∂+=00 (4) or equivalently j ij i u A u ~~=, where⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=100232221131211a a a a a a A ij (5)is an affine transformation obtained by linearising the homography about 0i u . This implies that each small image patch undergoes an affine transformation, and justifies the use of SIFT features which are partially invariant under affine change.Once features have been extracted from all n images (linear time), they must be matched. Since multiple images may overlap a single ray, each feature is matched to its k nearest neighbours in feature space (we use k= 4). This can be done in O （n n log ） time by using a k-d tree to find approximate nearest neighbours. A k-d tree is an axis aligned binary space partition, which recursively partitions the feature space at the mean in the dimension with highest variance.3. Image MatchingAt this stage the objective is to find all matching (i.e. overlapping) images. Connected sets of image matches will later become panoramas. Since each image could potentially match every other one, this problem appears at first to be quadratic in the number of images. However, it is only necessary to match each image to a small number of overlapping images in order to get a good solution for the image geometry.From the feature matching step, we have identified images that have a large number of matches between them. We consider a constant number m images, that have the greatest number of feature matches to the current image, as potential image matches (we use m= 6). First, we use RANSAC to select a set of inliers that are compatible with a homographybetween the images. Next we apply a probabilistic model to verify the match.3.1 Robust Homography Estimation using RANSACRANSAC (random sample consensus) is a robust estimation procedure that uses a minimal set of randomly sampled correspondences to estimate image transformation parameters, and finds a solution that has the best consensus with the data. In the case of panoramas we select sets of r= 4 feature correspondences and compute the homography H between them using the direct linear transformation(DLT) method. We repeat this with n= 500 trials and select the solution that has the maximum number of inliers (whose projections are consistent with H within a tolerance ∈pixels). Given the probability that a feature match is correct between a pair of matching images (the inlier probability) is i p , the probability of finding the correct transformation after n trials isH p (is correct ）n r i p ))(1(1--= (6)After a large number of trials the probability of finding the correct homography is very high. For example, for an inlier probability i p = 0.5, the probability that the correct homography is not found after 500 trials is approximately 14101-⨯.RANSAC is essentially a sampling approach to estimating H. If instead of maximising the number of inliers one maximises the sum of the log likelihoods, the result is maximum likelihood estimation (MLE). Furthermore, if priors on the transformation parameters are available, one can compute a maximum a posteriori estimate (MAP). These algorithms are known as MLESAC and MAPSAC respectively.3.2 Probabilistic Model for Image Match VerificationFor each pair of potentially matching images we have a set of feature matches that are geometrically consistent (RANSAC inliers) and a set of features that are inside the area of overlap but not consistent (RANSAC outliers). The idea of our verification model is to compare the probabilities that this set of inliers/outliers was generated by a correct image match or by a false image match.For a given image we denote the total number of features in the area of overlap f n andthe number of inliers i n . The event that this image matches correctly/incorrectly is represented by the binary variable {}1,0∈m . The event that the th i feature match (){}1,0∈i f is an inlier/outlier is assumed to be independent Bernoulli, so that the total number of inliers is Binomial()()()1:1,;1p n n B m f p f i n f == (7) ()()()0:1,;0p n n B m f p f i n f == (8) where 1p is the probability a feature is an inlier given a correct image match, and 0p is the probability a feature is an inlier given a false image match. The set of feature match variables (){}f i n i f ,...,2,1,= is denoted ()f n f:1 . The number of inliers ()∑==f ni i i f n 1 and B(.) is the Binomial distribution()()x n p p x n x n p n x B x ---=1!!!),;( (9) We choose values 6.01=p and 1.00=p . We can now evaluate the posterior probability that an image match is correct using Bayes’ Rule()()()()()()()f f f n n n f p m p m f p f m p :1:1:1111==== (10) ()()())1(1)0(011):1(:1====+=m p m f p m p m fp f f n n (11) We accept an image match if min ):1()1(p f m p f n >= ()()()1110,;1);(min 01,->==<p m p p n n B m p p n n B accept f i f i reject(12)Choosing values 610)1(-==m p and 999.0min =p gives the conditionf i n n βα+> (13)for a correct image match, where α= 8.0 and β = 0.3.Though in practice we have chosen values for 0p ,1p ,()0=m p ,()1=m p , and min p , they could in principle be learnt from the data. For example, 1p could be estimated by computing the fraction of matches consistent with correct homographies over a large dataset.Once pairwise matches have been established between images, we can find panoramic sequences as connected sets of matching images. This allows us to recognise multiple panoramas in a set of images, and reject noise images which match to no other images (see figure (2)).(a) Image 1 (b) Image 2(c) SIFT matches 1 (d) SIFT matches 2(e) RANSAC inliers 1 (f) RANSAC inliers 2(g) Images aligned according to a homographyFigure 1. SIFT features are extracted from all of the images. After matching all of the features using a k-d tree, the m images with the greatest number of feature matches to a given image are checked for an image match. First RANSAC is performed to compute the homography, then a probabilistic model is invoked to verify the image match based on the number of inliers. In this example the input images are 517×374 pixels and there are 247 correct feature matches.(a) Image matches(b) Connected components of image matches(c) Output panoramasFigure 2. Recognising panoramas. Given a noisy set of feature matches, we use RANSAC and a probabilistic verification procedure to find consistent image matches (a). Each arrow between a pair of images indicates that a consistent set of feature matches was found between that pair. Connected components of image matches are detected (b) and stitched into panoramas (c). Note that the algorithm is insensitive to noise images that do not belong to a panorama (connected components of size 1 image).4. Bundle AdjustmentGiven a set of geometrically consistent matches between the images, we use bundle adjustment to solve for all of the camera parameters jointly. This is an essential step as concatenation of pairwise homographies would cause accumulated errors and disregard multiple constraints between images, e.g., that the ends of a panorama should join up. Images are added to the bundle adjuster one by one, with the best matching image (maximum number of consistent matches) being added at each step. The new image is initialised with the same rotation and focal length as the image to which it best matches. Then the parameters are updated using Levenberg-Marquardt.The objective function we use is a robustified sum squared projection error. That is, each feature is projected into all the images in which it matches, and the sum of squared image distances is minimised with respect to the camera parameters. Given a correspondence lj k i u u (k i u denotes the position of the kth feature in image i ), the residual isk ij k i k ij p u r -= (14)where k ij p is the projection from image j to image i of the point corresponding to k i ul j j T j i i k ij u K R R K p ~1~-= (15)The error function is the sum over all images of the robustified residual errors()()()∑∑∑=∈∈=ni i j j i f k k ijr h e 1,ι (16) Where n is the number of images, I(i ) is the set of images matching to image i , ),(j i f is the set of feature matches between images i and j . We use a Huber robust error function()⎪⎩⎪⎨⎧≥-<=σσσσx if x x if x x h 222 (17) This error function combines the fast convergence properties of an 2L norm optimisation scheme for inliers (distance less than σ), with the robustness of an 1L norm scheme for outliers (distance greater than σ). We use an outlier distance ∞=σ during initialisation and 2=σ pixels for the final solution.This is a non-linear least squares problem which we solve using the Levenberg-Marquardt algorithm. Each iteration step is of the formr J C J J T p T 11)(--+=Φλ (18)Where Φ are all the parameters, r the residuals and Φ∂∂=r J . We encode our prior belief about the parameter changes in the (diagonal) covariance matrix p C⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡=2222200000000000000000000θθθθσσσσσf p C (19) This is set such that the standard deviation of angles is 16πθσ= and focal lengths10-=f f σ(where -f is the mean of the focal lengths estimated so far). This helps in choosing suitable step sizes, and hence speeding up convergence. For example, if a spherical covariance matrix were used, a change of 1 radian in rotation would be equally penalised as a change of 1 pixel in the focal length parameter. Finally, the λ parameter is varied at each iteration to ensure that the objective function of equation 16 does in fact decrease. The derivatives are computed analytically via the chain rule, for example1~~1i k ij k ij k ij i kij p p p p θθ∂∂∂∂=∂∂ (20) where [][]⎥⎦⎤⎢⎣⎡--=∂∂=∂∂22~1001z y z z x z z y x z y z x p p k ijk ij (21) and l j j j i i i i k iju K R R K p ~111~-∂∂=∂∂θθ (22) [][]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=∂∂=∂∂⨯⨯010********i i e e R i i i θθθθ (23) 4.1 Fast Solution by Direct Computation of the Linear SystemSince the matrix J is sparse, forming J J T by explicitly multiplying J by its transpose is inefficient. In fact, this would be the most expensive step in bundle adjustment, costing O(2MN ) for an N M ⨯ matrix J (M is twice the number of measurements and N is the number of parameters). The sparseness arises because each image typically only matches to a small subset of the other images. This means that in practice each element of J J T can be computed in much fewer than M multiplications()()1,-Φ∈=Φ∂∂Φ∂∂=∑C r r J J j k ij T j i f k i kij ij T(24)i.e., the inverse covariance between cameras i and j depends only on the residuals of feature matches between i and j .Similarly, r J T need not be computed explicitly, but can be computed via()()()k ij T n i i I j j i f k i k ij iT r r r J ∑∑∑=∈∈Φ∂∂=1, (25)In both cases each summation would require M multiplications if each feature matched to every single image, but in practice the number of feature matches for a given image is much less than this. Hence each iteration of bundle adjustment is O ()3N , which is the cost of solving the N N ⨯linear system. The number of parameters N is 4 times the number of images, and typically M is around 100 times larger than N .5. Automatic Panorama StraighteningImage registration using the steps of sections 2 - 4 gives the relative rotations between the cameras, but there remains an unknown 3D rotation to a chosen world coordinate frame. If we simply assume that R=I for one of the images, we typically find a wavy effect in the output panorama. This is because the real camera was unlikely to be perfectly level and un-tilted. We can correct this wavy output and automatically straighten the panorama by making use of a heuristic about the way people typically shoot panoramic images. The idea is that it is rare for people to twist the camera relative to the horizon, so the camera X vectors (horizontal axis) typically lie in a plane (see figure 4). By finding the null vector of the covariance matrix of the camera X vectors, we can find the “up -vector” u(normal to the plane containing the camera centre and the horizon)00=⎪⎭⎫ ⎝⎛∑=u X X T i n i i (26) Applying a global rotation such that up-vector u is vertical (in the rendering frame) effectively removes the wavy effect from output panoramas as shown in figure 4.Figure 3. Finding the up-vector u. A good heuristic to align wavy panoramas is to note that people rarely twist the camera relative to the horizon. Hence despite tilt (b) and rotation (c), the camera X vectors typically lie in a plane. The up-vector u(opposite to the direction of gravity) is the vector normal to this plane.(a) Without automatic straightening(b) With automatic straighteningFigure 4. Automatic panorama straightening. Using the heuristic that users rarely twist the camera relative to the horizon allows us to straighten wavy panoramas by computing the up-vector (perpendicular to the plane containing the horizon and the camera centre).6. Gain CompensationIn previous sections, we described a method for computing the geometric parameters (orientation and focal length) of each camera. In this section, we show how to solve for a photometric parameter, namely the overall gain between images. This is set up in a similar manner, with an error function defined over all images. The error function is the sum of gain normalised intensity errors for all overlapping pixels()()()()211,~~21∑∑∑===∈-=n i nj u H u j i R u j j j i i i j IJ i i u I g u I g e (27) where i g ,j g are the gains, and ()j i R , is the region of overlap between images i and j .In practice we approximate ()i u I by the mean in each overlapping region ij I -()()()∑∑∈∈-=j i R u j i R u i i ij i i u I I ,,1 (28)This simplifies the computation and gives some robustness to outliers, which might arise due to small misregistrations between the images. Also, since 0=g is an optimal solution to the problem, we add a prior term to keep the gains close to unity. Hence the error function becomes()⎪⎪⎭⎫ ⎝⎛-+⎪⎭⎫ ⎝⎛-=--==∑∑2222111/21g i N ji j ij i n i n j ij g I g I g N e σσ (29) where ()j i R N ij ,=equals the number of pixels in image i that overlap in image j . The parameters N σ and g σare the standard deviations of the normalised intensity error and gain respectively. We choose values {}()25500.10 ∈=I N σ and 1.0=g σ. This is a quadratic objective function in the gain parameters g which can be solved in closed form by setting the derivative to 0 (see figure 5).(a) Half of the images registered(b) Without gain compensation(c) With gain compensation(d) With gain compensation and multi-band blendingFigure 5. Gain compensation. Note that large changes in brightness between the images are visible if gain compensation is not applied (a)-(b). After gain compensation, some image edges are still visible due to unmodelled effects such as vignetting (c). These can be effectively smoothed out using multi-band blending(d).7. Multi-Band BlendingIdeally each sample (pixel) along a ray would have the same intensity in every image that it intersects, but in reality this is not the case. Even after gain compensation some image edges are still visible due to a number of unmodelled effects, such as vignetting (intensity decreases towards the edge of the image), parallax effects due to unwanted motion of the optical centre, misregistration errors due to mismodelling of the camera, radial distortion andso on. Because of this a good blending strategy is important.From the previous steps we have n images (){}()n i y x I i 1,∈ which, given the known registration, may be expressed in a common (spherical) coordinate system as ()φθ,i I . In order to combine information from multiple images we assign a weight function to each image ()()()y x y x W ωω=, where ()x ω varies linearly from 1 at the centre of the image to 0 at the edge. The weight functions are also resampled in spherical coordinates ()φθ,i W . A simple approach to blending is to perform a weighted sum of the image intensities along each ray using these weight functions()()()()∑∑===ni i ni i i linear W W I I 11,,,,φθφθφθφθ (30)Where ()φθ,linear I is a composite spherical image formed using linear blending. However, this approach can cause blurring of high frequency detail if there are small registration errors (see figure 7). To prevent this we use the multi-band blending algorithm of Burt and Adelson . The idea behind multi-band blending is to blend low frequencies over a large spatial range, and high frequencies over a short range.We initialise blending weights for each image by finding the set of points for which image i is most responsible()()()⎩⎨⎧==others W ifW W j ji i 0,max arg ,1,maxφθφθφθ (31) i.e.()φθ,max i W is 1 for ()φθ, values where image i has maximum weight, and 0 where some other image has a higher weight. These max-weight maps are successively blurred to form the blending weights for each band.A high pass version of the rendered image is formed()()()φθφθφθσσ,,,i i i I I B -= (32) ()()()φθφθφθσσ,,,g I I i i *= (33)Where ()φθσ,g is a Gaussian of standard deviation σ, and the * operator denotes convolution. ()φθσ,B represents spatial frequencies in the range of wavelengths []σλ,0∈. We blend this band between images using a blending weight formed by blurring the max-weight map for this image()()()φθφθφθσσ,,,maxg W W i i *= (34) where ()φθσ,i W is the blend weight for the wavelength []σλ,0∈ band. Subsequent frequency bands are blended using lower frequency bandpass images and further blurring the blend weights, i.e. for 1≥k()()i k i k ik I I B σσσ11++-= (35)()σσσ'+*=g I I i k ik 1 (36)()σσσ'*=+g W W i k ik 1 (37)where the standard deviation of the Gaussian blurring kernel ()σσ12+='k is set such thatsubsequent bands have the same range of wavelengths.For each band, overlapping images are linearly combined using the corresponding blend weights ()()()()∑∑===ni i k ni i k i k multik W W B I 11,,,,φθφθφθφθσσσσ (38)This causes high frequency bands (small σk ) to be blended over short ranges whilst low frequency bands (large σk ) are blended over larger ranges (see figure (6)).Note that we have chosen to render the panorama in spherical coordinates θ, φ.In principle one could choose any 2-dimensional parameterisation of a surface around the viewpoint for rendering. One good choice would be to render to a triangulated sphere, constructing the blending weights in the image plane. This would have the advantage of uniform treatment of all images, and it would also allow easy resampling to other surfaces (in graphics hardware). Note that the θ, φ parameterisation suffers from singularities at the poles.Algorithm: Automatic Panorama StitchingInput : n unordered imagesI. Extract SIFT features from all n imagesII. Find k nearest-neighbours for each feature using a k-d treeIII. For each image:(i) Select m candidate matching images that havethe most feature matches to this image(ii) Find geometrically consistent feature matchesusing RANSAC to solve for the homographybetween pairs of images(iii) Verify image matches using a probabilistic modelIV . Find connected components of image matchesV . For each connected component:(i) Perform bundle adjustment to solve for therotation 1θ,2θ,3θ and focal length f of all cameras(ii) Render panorama using multi-band blendingOutput : Panoramic image(s)(a) Original images and blended result(b) Band 1 (scale 0 to σ)(c) Band 2 (scaleσto2σ)(d) Band 3 (scale lower than2σ)B(θ, φ)for k= 1,2,3 are shown on the left, with the Figure 6. Multi-band blending. Bandpass images σkW(θ, φ) shown on the right. Initial blending weights are assigned to 1 corresponding blending weights σkwhere each image has maximum weight. To obtain each blending function, the weights are blurred at spatial frequency σ and bandpass images of the same spatial frequency are formed. The bandpass images are blended together using weighted sums based on the blending weights (Note: the blending widths have been exaggerated for clarity in these figures).(a) Linear blending (b) Multi-band blendingFigure 7. Comparison of linear and multi-band blending. The image on the right was blended using multi-band blending using 5 bands and σ= 5 pixels. The image on the left was linearly blended. In this case matches on the moving person have caused small misregistrations between the images, which cause blurring in the linearly blended result, but the multi-band blended image is clear.(a) (b)Figure 8. Stitching with rotation and zoom. Our use of invariant features make stitching possible despite rotation, zoom and illumination changes in the input images. Here the inset images at the base and tip of the tower are 4 times the scale of the other images.(a) Number of feature matches vs κ (b) Example of stitching 2 (c) Detail from the images with κ=−0.5(the previous figure sho-full test set contains 44 wing misregistrationimages of this scene)Figure 9. Stitching with radial distortion. This figure shows the effects on stitching of first order radial distortion ()x x k x 21+=' with κ in the range []5.0,5.0-∈k (the image height is normalised to unit length). Note that radial distortion is not modelled in our algorithm. We used a test sequence of 44 images and applied radial distortion with 20 values of κ. Examples of the distorted images are given in figures (d)-(h). To evaluate the performance of stitching we counted the number of consistent matches after RANSAC, the results are shown in figure (a). Although the number of matches per feature dropped by around a third in the worst case, the number of correct feature matches was still high(around 500 per image), so the images could still be successfully matched. Nevertheless radial distortion causes visible artifacts in rendering as shown in figures (b)-(c), and correcting for this in the bundle adjustment and rendering stages would be important for high quality panorama stitching.8. ResultsFigure 2 shows typical operation of the panoramic recognition algorithm. A set of images containing 4 panoramas and 4 noise images was input. The algorithm detected connected。