InternationalJournalofComputerVision64(2/3),125–141,2005 (cid:1)c 2005SpringerScience+BusinessMedia,Inc..ManufacturedinTheNetherlands. Scale-Space Image Analysis Based on Hermite Polynomials Theory SHERIFMAKRAM-EBEID* ANDBENOITMORY PhilipsMedicalSystemsResearchinParis, 51rueCarnot,B.P.301, F-92156,SURESNESCedex,France sherif.makram-ebeid@philips.com benoit.mory@philips.com ReceivedOctober8,2003;RevisedNovember11,2004;AcceptedNovember11,2004 FirstonlineversionpublishedinJune,2005 Abstract TheHermitetransformallowstolocallyapproximateanimagebyalinearcombinationofpolynomials. For a given scale σ and position ξ, the polynomial coefficients are closely related to the differential jet (set of partial derivatives of the blurred image) for the same scale and position. By making use of a classical formula due to Mehler (late 19th century), we establish a linear relationship linking the differential jets at two different scalesσ andpositionsξ involvingHermitepolynomials.Formulti-dimensionalimages,anisotropicexcursionsin scale-spacecanbehandledinthisway.Patternregistrationandmatchingapplicationsaresuggested. WeintroduceaGaussianwindowedcorrelationfunctionK(v)forlocallymatchingtwoimages.Whentakingthe mutual translation parameter v as an independent variable, we express the Hermite coefficients of K(v) in terms of the Hermite coefficients of the two images being matched. This new result bears similarity with the Wiener- KhinchintheoremwhichlinkstheFouriertransformoftheconventional(flat-windowed)correlationfunctionwith theFourierspectraoftheimagesbeingcorrelated.Comparedtotheconventionalcorrelationfunction,oursismore suitedformatchinglocalizedimagefeatures. Numerical simulations using 2D test images illustrate the potentials of our proposals for signal and image matchingintermsofaccuracyandalgorithmiccomplexity. Keywords: Hermitetransform,anisotropicscale-space,patternmatching,affineregistration,Gaussianwindowed correlation,localcorrelation 1. Introduction ilyofwavelettransformsaimatfulfillingsimilarneeds (Daubechies,1992;Mallat,1998). One of the key issues in signal and image process- In the present article, we study orthogonal polyno- ing is to identify appropriate data representations to mialrepresentationsthataresuitedtorepresentsignal suit a given range of applications. The widest spread andimagedatawithinaGaussianwindow.Suchawin- such representation is the Fourier transform which dowisuniqueinthesensethatitisoptimalwithrespect is quite suited, for example, to capture the oscilla- to scale-space axioms and is dimensionally separable torybehaviourofsignalsandnoise(Bracewell,1999; (Duitsetal.,2002;Floracketal.,1996).Furthermore, Schro¨der and Blume, 2000). When local features or asengineersandphysicistsknow,itoffersthebestsi- patterns have to be analysed or processed, Gaussian multaneousfrequency-spacelocalizationcompromise windowed Fourier or Gabor transforms are preferred (Arfken,1985). (KrugerandSommer,2002;Lee,1996).Alargefam- For a fixed Gaussian window, the corresponding Hermitepolynomialrepresentationisstronglyrelated *Towhomcorrespondenceshouldbeaddressed. to the biologically motivated Gabor representation 126 Makram-EbeidandMory (Lee, 1996). It has advantages over Gabor wavelets and because it is more convenient for handling rotations, (cid:1) translationsandscalechangesaswillbeshowninthe +∞ c = w(x)(p (x))2dx. (2) next sections. As with Gabor or with other wavelet i i −∞ transforms,systematicanalysisorprocessingofacom- plete signal or image requires an array of overlap- pingwindows(KrugerandSommer,2002;Lee,1996; Theweightfunctionw(x)definesan L metricinthe 2 Martens, 1997). In general, this should be preferably vector-space of functions f(x) for which any distinct doneforarangeofdifferentscales(i.e.differentsizes pairofpolynomialsinthebasisareorthogonal.Adirect of Gaussian windows in our case) within a multi- consequenceofthisisthatonecandefineascalarprod- resolution pyramid. A very powerful tool for dealing uctbetweenanytwofunctions f(x)andg(x)through withsuchwindowarraysisthewaveletframestheory Parseval’stheorem(Arfken,1985) (Daubechies,1992;Lee,1996;Mallat,1998;Martens, 1997). In this article, we concentrate on the local (cid:1) +∞ (cid:3)∞ f g aspects that constitute essential building blocks of (cid:3)f,g(cid:4) ≡ w(x)f(x)g(x)dx = i i w a multiple window multi-resolution pyramid that we −∞ i=0 ci wishtodevelopinthefuture. where 1.1. OrthogonalPolynomialsforSignal (cid:1) +∞ andImageProcessing g = w(x)g(x)p (x)dx. (3) i i −∞ We consider signals and images as scalar valued functions of one or several variables. We start here The weighted scalar product (cid:3)f,g(cid:4) can therefore be w by dealing with one dimensional signals or, equiv- expressed either as an integral in x-space or else in alently, functions of one real variable. Orthogonal terms of polynomial transform coefficients f and g . i i polynomials provide a convenient tool for locally ap- ThediscretesumintherightmostmemberofEq.(3)de- proximating such functions (Dunkl and Xu, 2001; fines an equivalent transform-space (index i-labelled) Szego¨,1967).Thus,oneapproximatesafunction f(x) scalarproduct;itscorrespondingmetricisreferredto, as a linear combination of polynomials p (x) having i inthisarticle,asthel -metric.Theinverseofthenor- 2 degreesigoingfrom0toinfinity.Theapproximation malization constants c introduce weights in the l - i 2 is performed by a least mean square technique that metric playing a role similar to that of w(x) in the consistsofminimizingtheweightedquadraticerrorQ x-space L -metriccorrespondingtothescalarproduct 2 definedby (cid:3)f,g(cid:4) . w (cid:2) (cid:4) Theaboveresultsallowtherepresentationofafunc- (cid:1) +∞ (cid:3)∞ 2 tion of a real variable by an infinite polynomial se- Q = w(x) f(x)− a p (x) dx (1) i i ries.Inpracticeonlyafinitenumberoftermsmaybe −∞ i=0 used.Skippingthehighranktermsgivesrisetotheso calledGibbsoscillationsphenomenon(Mallat,1998). where ai are the polynomial coefficients and where For the particular case of Hermite polynomials dealt w(x) is a non-negative weight function. Regions of withinthisarticle,largepolynomialorderasymptotic x for which w(x) is large are those for which one expansions can be found in Szego¨ (1967). Those can wishes the approximation to be accurate even when beusedtoshowthatmissinghighranktermshavean keeping a small number of terms. The coefficients osc√illatory behaviour with amplitudes behaving like ai arethosewhichminimizethequadraticerrorQ.For 1/ w(x). In other words, the square of the approx- apolynomialbasiswhichisorthogonalfortheweight imation error is inversely proportional to the weight functionw(x),theyaregivenby functionw(x).Thisiswhatonewouldintuitivelyex- (cid:1) pectfromEq.(1);i.e.thesquareoftheapproximation +∞ errorislargewhenitsconstrainingweightfactorw(x) a = f /c where f = w(x)f(x)p (x)dx i i i i i issmall. −∞ Scale-SpaceImageAnalysisbasedonHermitePolynomialsTheory 127 1.2. HermiteTransformandRelatedScale-Space transformisdirectlyrelatedtothederivativesof F(x) DifferentialStructure whichdefinethedifferentialstructureoffunction f(x) atscaleσ. Asmentionedearlier,theGaussianwindowisanatural choice for image processing. For one dimension, we useforw(x)anormalizedwindowofarbitrarysizeσ 2. BasicPropertiesofHermiteTransform andcentroidξ (cid:5) (cid:6) In this main section, we study the basic properties of 1 (x −ξ)2 theHermitetransformthatcanbeusefulforsignaland w(σ,ξ;x)= √ exp − . (4) 2πσ2 2σ2 imageprocessing.Wemakeuseofthetheoryoforthog- onal polynomials which is well established (Szego¨, 1967)andwhichhasbeenextendedtoanynumberof Thecorrespondingorthogonalbasisisthatoft√heHer- dimensions(DunklandXu,2001). mitepolynomials(Szego¨,1967) H((x−ξ)/σ 2)for i which the normalization constants c defined in Eq. i (2) is given by c = 2ii!. The Hermite transform of 2.1. DimensionalSeparability i anarbitraryfunction f(x)isgivenbythecoefficients anditsComputationalAdvantages f defined in Eq. (2) with the above special weight i function and with pi(x) replaced by the above Her- Ad-dimensionalimageisdefinedhereasareal-valued mitepolynomials.Hermitecoefficients fi canalsobe function f(x) of the d-dimensional position vector mathematically written in terms of derivatives of the x = (x ,x ,...,x )T. As for the 1D case, one can 1 2 d function F = w ⊗ f obtained by convolving func- locallyanalyse f(x)withinanisotropicGaussianwin- tion f(x) with a normalized Gaussian kernel of size dow of size σ (representing scale) around any win- σ.ByusingthebasicpropertiesofHermitepolynomi- dow centre ξ = (ξ ,ξ ,...,ξ )T. The isotropic d- 1 2 d als(Szego¨,1967),itisasimpleexercisetoshowthat dimensional window can be considered as a product (KoenderinkandVanDoorn,1992) ofone-dimensionalwindows.Likewise,theorthogonal (cid:7) polynomialbasisismadeupofproductsof1DHermite fi(σ,ξ)= 2i/2σiddixFi (cid:7)(cid:7)(cid:7) . (5) p(DoluynnkolmaniadlsX,wu,h2ic0h01a)r.eEsahcohwonnetooffotrhmosaecpoomlypnloemteiaselst x=ξ is identified by a multi-index or d-dimensional vec- torindex I = (i ,i ,...,i )whereallcomponentin- A particular class of functions f is that of the poly- 1 2 d dicesi arenon-negativeintegers.Thecorresponding nomials. In this case, the blurred function F is also a k d-dimensional Hermite transform coefficient f can polynomial and the resulting Hermite expansion will I nowbewrittenas includeafinitenumberofnon-vanishingterms.There- sultsobtainedinthisarticlecanthenbeconsideredas 1 exacttoolstomanipulatepolynomialapproximations. f (σ,ξ)= I (2πσ2)d/2 Asimilarremarkextendstothemulti-dimensionalcase (cid:1) (cid:2) (cid:4) (cid:9)x −ξ(cid:9)2 (cid:8) (cid:9) (seeSection3below). exp − P x f(x)dx √Iftheindependentvariablexisexpressedinunitsof x∈Rd 2σ2 I σ 2,thepre-factor2i/2σi inEq.(5)canbeomitted. Equation(5)willbeusedlatertoderiveotherresults. where Itcanbeused,inpractice,toobtainlow-orderHermite (cid:5) (cid:6) bcoluerfrfiecdievnetrssifoinwFit(hx)i o pand A(p) with j =0arereplacedbyzero. i,j−1 F(σ 1−t2;ξ +tv) For image dimensions d > 2 any rotation can be (cid:5) (cid:6) decomposedintoasuccessionofin-planerotationsin- =(cid:3)∞ (σt)i H √v F(i)(σ;ξ) volving two coordinates at a time (Gallier, 2000). A 2i/2i! i σ 2 i=0 straightforward extension of the above approach can (cid:5) (cid:6) (cid:3)∞ ti f (σ;ξ) v thereforedealwithmulti-dimensionalrotationsinHer- = i H √ (17) mitedomain. 2ii! i σ 2 i=0 130 Makram-EbeidandMory where F(γ;η) stands for the signal f Gaussian If one is interested in polynomial approximations smoothed with kernel σ = γ and evaluated at x = η not exceeding a maximum degree N, all coefficients and F(i)(γ;η) is the corresponding ith derivative. By of order larger than N can be replaced by zeros and takingtclosetounity,onereconstructsthesignalwith- the above relations are expressed as a linear relation outblurring. betweenthe(N +1)-dimensionalHermitecoefficient vectors for the two windows. The matrix incurred in Mehler’s Formula Related Signal Processing this transformation is an upper triangular one where LEMMA. Equation 17 implies that, given a factor any row j is obtained from the first one (j = 0) by a t with|t|<1,iftheHermitecoefficientsaremodified shift and a scalar multiplication by αj. In the special bymultiplyingordericoefficientbyti;thentheresult- case where f(x) is a polynomial of degree N, the ingHermiteseriesexpansionyieldsavers√ionofsignal resultingformulaeareexact. f blurredwithaGaussiankernelofsizeσ 1−t2and When the function f(x) is not a pure polynomial, thenzoomedbyafactor1/t aboutthewindowcentre. one still wishes to use Eq. (18) with all terms of or- AnequivalentformulationwasderivedbyMartens derlargerthanNneglected.Thisisjustifiedaslongas (1992) for deblurring isotropic Gaussian blur in im- theapproximationerrorsduetohigherordertermsare ages. withinacceptablelimits.Figure1illustratesthevalid- itycondition.Thecoefficients fj+i(σ;ξ)correspondto New Analytical Expression for 1D Scale-Space theGaussianwindowW1(x)withkernelsizeσandcen- Transformations. Florack et al. (1996) have anal- treatx =ξwhereascoefficients fj(σα;ξ+βv)corre- ysedtheinfluenceofsimultaneoustranslationandscale spondtotheGaussianwindowW2(x)withkernelsize reduction excursions on scale-space local jets. They σαandcentreatx =ξ+βv.Asdiscussedattheendof havestudiedtheconstraintsonsuchexcursionsthaten- Section 1.1, approximation errors within window W1 sureconvergenceoftheirTaylorexpansions.Mehler’s canbeexpectedtobesmallfortherangeofxvaluesfor formula allows us to more specifically put the corre- whichW1(x)islargerthanapositivefraction(cid:9).W1(ξ) spondingresultsinanexplicitandcompactanalytical ofitspeakvalueW1(ξ)(with0<(cid:9) <1,e.g.(cid:9) =0.1). form provided appropriate pairs of “conjugate zoom- Likewise, the validity range within window W2 can ing and blurring operators” are applied. To do so, we be defined by the condition W2(x) > (cid:9).W2(ξ +βv). recall that Eq. (5) allows to mathematically express Since Eq. (18) expresses Hermite coefficients for W2 HermitecoefficientsintermsofderivativesofaGaus- sianblurredfunction.WeapplyEq.(17)withsignalf replaced by its jth derivative and the values of F and ofitsderivativesexpressedintermsofHermitecoeffi- cientstoget f (σα;ξ +βv) j (cid:2) (cid:5) (cid:6) (cid:4) (cid:3)∞ βi v =αj 2ii!Hi σ√2 fj+i(σ;ξ) (18) i=0 √ where α = 1−t2 and β = t so that α2 +β2 = 1. ThisallowstoexpressHermitecoefficientsforaGaus- Figure 1. Illustrating the validity condition for the 1D scale- sianwindowtranslatedtox =ξ+βvandreferredtoa space transformation deduced from Mehler’s formula. Given windowofsizeσαintermsoftheHermitecoefficients the Hermite coefficients fi(σ;ξ) corresponding to the Gaussian for an original window centred at x = ξ and of size cwoienfdfiocwienWts1(fxj)(,σEαq;.ξ(1+8)βavl)locwosrrethsepocnodminpguttaotiothneoGfatuhsesiHanerwmiinte- σ.ItisusefultonotethattheGaussiankernelsizesσ dow W2(x). For W2(x), the kernel size σα must be smaller than andσαarescales(i.e.sticklengths)inunitsofwhich that the kernel size σ of W1(x) (i.e. α < 1). The relative trans- lationβv aswellastherelativescalefactorα mustbesuchthat the independent variable x may be measured. Since α ≤1,theHermiteexpansionfortheGaussiankernel a{xsm|Wall1(pxo)si>tiv(cid:9)e.Wfra1c(tξi)o}n⊂of{uxn|ityW(2(cid:9)(x=)>0.1(cid:9).sWay2)(.ξI+noβthve)}rwwohredres,(cid:9)thies σα maybeconsideredasazoomedversionofthatof acceptableHermiteexpansionaccuracyrangeforW2(x)shouldbe theGaussiankernelσ. includedinthatofW1(x) Scale-SpaceImageAnalysisbasedonHermitePolynomialsTheory 131 √ √ intermsofthoseofW1,itsvalidityrequiresthattheset size σα1 1−t2 and σα2 1−t2 respectively. Par- {x | W2(x) > (cid:9).W2(ξ +βv)} be included within the seval theorem allows writing the L2 (x-space) dif- set{x | W1(x) > (cid:9).W1(ξ)}.Inotherwords,thegood ference norm of the difference of blurred functions accuracyrangeofW2shouldbeincludedinthatofW1. in terms of a modified index-labelled l2 difference Itisclear,inparticular,thatwindow W2 couldnotbe norm: a pure translated version of W as this would always 2 (cid:1) vchiooliacteeothfe(cid:9)ab>ov0e.aTchciusriascywrhayngtreainnsclalutisoionnmruulsetfaolwraanyys √1π +∞ (cid:8)e−u22F(ξ +β1v1+α1tu) u=−∞ be accompanied by zooming (i.e. reducing the kernel (cid:9) size).Thesameremarksalsoexplainwhyonecannot −e−u22G(ξ +β2v2+α2tu) 2du deaAlwuistehfsuilgnsapleschiarlinkcaasgee(iis.e.thtaaktinogfα1D>1h)o.methety ∼=(cid:3)N t2i (cid:8)f (σα ;ξ +β v ) 2ii! i 1 1 1 (Gallier, 2000) which is an affine transformation for i=0 (cid:9) whichthetranslationparameterviszeroandforwhich −g (σα ;ξ +β v ) 2. (22) i 2 2 2 thecentreoftheGaussianwindowisafixedpoint.All oddorderHermitepolynomialsinEq.(18)vanishand, makinguseoftheidentityβ2 =1−α2,oneisleftwith thesimplerexpression 2.4. NewHermite-Domain1DWindowed (cid:2) (cid:4) (cid:3)∞ (α2−1)i Cross-Correlation fj(σα;ξ)=αj 22ii! fj+2i(σ;ξ) . (19) i=0 In the previous section, we have proposed a method to evaluate the mean square difference of two signals formatchingpurposes.Wenowinvestigatethecorre- New 1D l -Matching based on a Further Use of 2 lation function for two signals versus relative trans- Mehler’s Formula. Take two signals f and g for lation and scale parameters. We attempt to do this in which Hermite coefficients are known in an original windowcentredatx =ξ andofsizeσ.Equation(18) such a way that the two signals being correlated play symmetricalroles.Furthermore,webuildthecorrela- canbeusedtocomparescaledandtranslatedversions tionfunctioninahierarchicalfashion(coarsetofine) of the two signals in Hermite coefficient domain. To usingthesmallestpossiblenumberofHermitecoeffi- do so, we operate a scale and translation change on f and g with parameters α = α ,β = β and v = v cientsofeachsignal.Tosimplifynotation,weassume i i i withi = 1 for f andi = 2 for g and with constraints that the Gaussian window is centred at the ordinates that α2 +β2 = 1. One can, furthermore, smooth and o√rigin and that the ordinates unit is chosen so that i i 2σ2 =1. zoom the two functions before comparing them. The Taketwo1Dfunctions f(x)andg(x)withHermite correspondingHermiteapproximationstoordernare transforms f andg .Ratherthanpredictinghoweach i i F(ξ +β v +α tu) ∼= (cid:3)N ti fi(σα1;ξ +β1v1) otifonthoersesctawloecchoaenffigcei,ewntessetutsdayreinmsteoaddifitehdebweihthavtriaonusrloaf- 1 1 1 i=0(cid:5) (cid:6)2ii! theexternalproducth(x,y)= f(x).g(y)ofthesetwo u functions. We thus generate a 2D function out of the × H √ (20) i σ 2 two1Dfunctions.The2DHermitecoefficientsofthis new2Dfunctionaregivenbyhi,j = figj.Letusnow and perform a (−θ) rotation from the (x,y) referential to therotated(u,v)sothat G(ξ +β v +α tu) 2 2 2 (cid:2) (cid:4) (cid:2) (cid:4)(cid:2) (cid:4) (cid:2) (cid:4) (cid:5) (cid:6) ∼=(cid:3)N tigi(σα2;ξ +β2v2)H √u (21) u = cosθ sinθ x andhence x 2ii! i σ 2 v −sinθ cosθ y y i=0 (cid:2) (cid:4)(cid:2) (cid:4) cosθ −sinθ u = . (23) where F and G are obtained from the functions f and sinθ cosθ v g by convolving with Gaussian low-pass kernels of 132 Makram-EbeidandMory Equation (12) tells us how to use the above hi,j co- u;theyareequivalenttoanaffinetransformationfrom efficients in the (x,y) referential in order to com- argumentxofftoargumentyofgwhichcanbewritten putethecoefficientsh(cid:11) inthe(u,v)referentialwith as 0,p the first index equal to zero (i = 0) and the sec- ond index equal to an arbitrary non-negative integer y =ζx +α−1v, with inverse x =ζ−1y−β−1v (j = p). From the remark in the last paragraph of Section 2.1, one sees that these h(cid:11) coefficients are (28) 0,p justtheHermitecoefficients of the1Dfunction K(v) obtained by Gaussian smoothing in the u direction, where ζ is a relative scaling parameter and α−1v (or i.e. −β−1v)asrelativetranslations. (cid:1) +∞ 1 K(v)= √ e−u2 f(α.u−β.v)g(β.u+α.v)du π AdvantagesofourHermite-DomainWindowedCor- −∞ (24) relation. Equations (25)–(28) provide new original features relative to earlier proposed correlation func- where α = cosθ and β = sinθ. In other words, the tions. They allow to deal with translation for any Hermite coefficients Kp of K(v) are given by h(cid:11)0,p givenscalechangewhencomputingtheGaussianwin- whichfromEq.(12)gives dowed cross-correlation (Eq. (24)). Furthermore, the Hermite transform of the cross-correlation function (cid:5) (cid:5) (cid:6) (cid:6) K =h(cid:11) =(cid:3)j=p (−1)j p sinjθcosp−jθ can be deduced explicitly from the Hermite coeffi- p 0,p j cients of the two functions f and g. When search- j=0 ing for a best match, our computational savings can · fj.gp−j. (25) be achieved by working in transform space and con- centratingonnoise-robustlow-ordercoefficients.Den Brinker (1993) have earlier reported related Laguerre Thefactthat K(v)isawindowedcross-correlationof transform domain correlation results. However, the the two functions f and g is particularly clear if one Ref.(DenBrinker,1993)LaguerrecounterpartofEq. choosesθ =π/4which,withachangeofindependent (25)isalessconvenientinfiniteseries(inplaceofthe variableintheintegrandofEq.(25),yields finite sum in Eq. (25)). In addition to its more com- (cid:5) (cid:6) (cid:13) (cid:1) (cid:14) (cid:15) pactform,ourHermitedomaincorrelationcanreadily τ 2 +∞ τ K √ = e−2u2 f u− extendtoseveraldimensionsbyvirtueofdimensional 2 π(cid:14)−∞ (cid:15) 2 separability. τ ×g u+ du. (26) 2 3. DealingwithanyNumberofImage The special case for θ = π/4 and with identical f Dimensions andgfunctionsdefinesa“windowedauto-correlation”. More generally, one may remark that the indepen- The dimensional separability of the Hermite trans- dent variable u in the integral of Eq. (25) is multi- form greatly facilitates dealing with any number of plied by a factor cosθ in the argument of f and by dimensions. A particularly useful feature is the rela- sinθ intheargumentofg.Theindependentvariablev tion with the differential structure of an image. The parameterises a translation of the argument of f rela- resultexpressedbyEq.(5)isreadilygeneralizedfora tive to that of g. For a fixed translation parameter, by d-dimensionalimagewithanisotropicGaussianwin- setting dow of size σ around a centre ξ. The Hermite coef- (cid:11) ficien(cid:16)t for a multi-index I = (i1,i2,...,id) of rank cosθ =α =1/(cid:11)1+ζ2 and n = dk=1ik isgivenby: sinθ =β =ζ/ 1+ζ2 (27) (cid:7) (cid:7) ∂nF (cid:7) the above rotation matrix equations relating (x,y) to fI(σ,ξ)=2n/2σn ∂xi1∂xi2...∂xid(cid:7)(cid:7) . (29) (u,v)canbereinterpretedbyeliminatingtheparameter 1 2 d x=ξ Scale-SpaceImageAnalysisbasedonHermitePolynomialsTheory 133 3.1. ReducingandExtendingtheNumberof comparingthemtogether.Todoso,welinearlyexpress TransformDimensions the pair of d-dimensional vectors (x,y)T in terms of a new pair of d-dimensional vectors (u,v)T. In this In Section 2.4, we have generated a 2D function newrepresentation,visafixedmutualtranslationvec- by taking the external product of two 1D functions. torparameterwhereasu isthenewindependentspace Thisprovedusefulforcomputingthelocalizedcross- variable. In order to apply the results of Sections 2.3 correlationanditsHermitetransform.Thisprocedure and2.4tothekthdimension,thecoordinatesx andy k k canbeextendedtoanynumberofdimensions.Takea ofthepositionalvectors x and y aretransformedinto function f(u)definedform-dimensionalvectorargu- u andv throughaunitarytransformationdefinedby k k mentsuandafunctiong(v)definedforn-dimensional (cid:2) (cid:4) (cid:2) (cid:4)(cid:2) (cid:4) positional vectors v, their product f(u).g(v) can be x α −β u treated as a function of the m +n-dimensional vec- yk = βk α k vk where αk2+βk2 =1 tor x = (u,v) = (u ,...,u ,v ,...,v ). It is not k k k k 1 m 1 n (30) hard to see that the Hermite coefficient of this product function for the m + n-dimensional multi-index I =(i1,i2,...,im+n)isequalto f(i1,...,im)g(im+1,...,im+n). aIfptahritsitiiosnd-monaetrifxorrealalltioimnabgeetwdeimenenthsieontwsok,(2w×e dge)-t The inverse situation is also of interest. In Section dimensionalvectors(x,y)T and(u,v)T intheform 2.1, for example, it was shown that f0,j are the 1D Hermite coefficients of f (y) obtained by Gaussian (cid:2) (cid:4) (cid:2) (cid:4)(cid:2) (cid:4) 0 smoothinga2Dfunction f(x,y)alongthex-axis.This x A −B u = (31) wasinstrumentalinderivingthenewHermitedomain y B A v correlationresultofEq.(25). whereAandBareddiagonalmatriceswithdiagonalel- 3.2. LocallyAffineWarpingisReadilyHandledin ementsgivenbyα andβ respectively.Recallingthat k k HermiteTransformDomainUsingourResults visthetranslationparameterandeliminatingthecom- monindependentspacevectorubyusingEq.(31),the TheresultsobtainedinSections2.3to2.4readilyex- correspondingimage-warpingtransformation(from x tendtoanynumberofdimensionsthroughthedimen- to y)isexpressedby sionalseparabilitypropertiesoftheHermitetransform. For any one of the image dimensions (labelled k), y = Zx + A−1v, withinversex = Z−1y−B−1v weapplyMehler’sformula(Eq.(17))tocomputeeach (32) oneofthemulti-dimensionalHermitecoefficients(Eq. (29))considered asfunctions of x .Thisallowsusto k linearlyexpressthedifferentialjetforgiventranslation where Z = BA−1 = A−1Bisadiagonald×d matrix and scale-reduction along x in terms of the original k withdiagonalelementsgivenbyζ = β /α andma- differential jet. This operation can be repeated for all k k k trix Z−1isitsinversewithdiagonalelementsgivenby image dimensions thus providing a generalization of 1/ζ =α /β . the results of Section 2.3. It is noteworthy that dif- k k k This, however, is not the most general image- ferent scale reductions (inverse of zooming factors) warping class that can be conveniently dealt with in canbeappliedtothedifferentdimensionsthusallow- Hermitetransformdomain.Asshownearlier,purero- ing anisotropic excursions in scale-space. As for the tationsareeasilyhandledinviewofthetensorproper- 1D case, neither pure translation nor size shrinkage tiesoftheHermitecoefficients.Onecanapplyrotation (equivalenttoscaleincrease)canbedirectlydealtwith operators R and R inx andin y respectivelyaround inHermitedomainasthiswouldviolatetheaccuracy 1 2 theirwindowcentressothatthegenericformoftrans- rangeinclusionprincipleexplainedinSection2.3and formation that can conveniently be handled is of the Fig.1.Inordertomatchimagesformoregeneraltypes form: of mutual warping laws using Hermite domain processing, different methods are possible. One method y = RTZR x + RTA−1v, withinverse consistsinapplyingdifferentzoomandtranslationto 2 1 2 eachofthed-dimensionalimages f(x)andg(y)before x = RTZ−1R y− RTB−1v. (33) 1 2 1 134 Makram-EbeidandMory If all axial scaling factors ζ are non-zero, the above 4. NumericalStudyforHermiteDomain k warping laws are invertible affine transformations RegistrationandMatching (Gallier,2000).Conversely,anyinvertibleaffinetrans- formationcanbeputinsuchformbyvirtueoftheSin- Wehaveperformednumericalsimulationstoillustrate gularValueDecompositionTheorem(Gallier,2000). theresultsreportedinthisarticleusing2Dimages.In Fig. 2(a), a 128×128 butterfly test-image is approx- imatedbyHermitepolynomialexpansionsuptorank N =64.TheHermitetransformimplementationused 3.3. Extendingourl NormandCorrelationResults forthispurposeisofthecontinuouswavelet-transform 2 toSeveralDimensions type (Daubechies, 1992; Mallat and Wavelet, 1998; Martens, 1997), and is based on the numerical evalu- GeneralizingEq.(17)toseveralvariablesisstraightfor- ation of the integrals involved in Eq. (8). This imple- ward.Itiseasilyshown,forexample,thatmultiplying mentation allows us to extract a bivariate polynomial eachHermitecoefficientbyafactortn where|t| < 1 approximation of the image valid for a region of the and where n is the coefficient’s rank (sum of the in- realx,yplanewithinaGaussiananalysiswindow(i.e. dicesi ),theresultingHermiteexpansionisaversion notonlyfordiscrete-valuedx,y coordinates).There- k of√f blurred with an isotropic Gaussian kernel of size sults of Sections 2 and 3 show how one can handle σ 1−t2 and zoomed about the window centre by a such polynomial approximants in an exact analytical factorof1/t.TheotherresultsofSections2.3and2.4 manner. In Figs. 2(a) and (e), the same butterfly im- extendtoseveraldimensionsbymakinguseofthedi- ages,havinggreyvaluesintheclosedinterval[0,255], mensionalseparabilityoftheHermitetransform.The isanalysedusingtwodifferentGaussiananalysiswin- genericclassofreferentialchangesofSection3.2can dowsschematizedbydottedcircleshavingtheircentres beeasilydealtwithinthismanner. at the Gaussian window centre and radii equal to 2σ Figure2. Thesamebutterfly128×128testimageisanalysedwithintwodifferentGaussiananalysiswindows(a)and(e)schematizedbythe dottedcircles.Bothcircleradiiareequalto2σwithσ =16.ThecorrespondingHermiteexpansionstomaximumtotalrankN =64areshown withdifferentblurandzoomparametertin(b),(c)and(d)startingfromthewindowin(a)andin(f),(g)and(h)startingfromthewindowin (e).Theabsolutegrey-leveldifference:between(b)and(f)isshownin(i),between(c)and(g)isshownin(j)andbetween(d)and(h)isshown in(k).

Scale-Space Image Analysis Based on Hermite Polynomials Theory

Similar Scale-Space Image Analysis Based on Hermite Polynomials Theory