Scale-Space Image Analysis Based on Hermite Polynomials Theory

17 Pages·2005·8.02 MB·English

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
[email protected]
[email protected]
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<f 4f(bxy).nHuomweervicear,lltyhidsenriuvminegria- PI(x)=(cid:10)d Hik x√k 2−σξ2k (6)
calmethodcannotbeextendedtocomputehigheror- k=1
der coefficients because of the occurrence of severe
roundingerrors.Computationofthesetofcoefficients and the corresponding polynomial approximation for
f , from Eq. (2), requires the knowledge of the non- f(x)inthewindowis
i
blurred function f(x) for x-values at which the win-
(cid:3)
dow function w(x) is not negligibly small. With this f(x)∼= f (σ,ξ)P (x)/c (7)
I I I
reservation in mind, Eq. (5) shows that the Hermite
I
128 Makram-EbeidandMory
where cI is the normalization constant which for a theorderof N2d whichismuchlarger.Thelargerthe
multi-index I =(i1,i2,...,id)isexpressedasaprod- imagedimension,themoresignificantisthesaving.
uctof1Dnormalizationfactors.ThesummationinEq. Another useful consequence of Eq. (8) is that f0,j
(7) extends over all multi-indices I one wishes to use can be interpreted as the 1D Hermite transform of
intheapproximation,allintegercomponentindicesik the function f0(y) obtained by Gaussian smoothing
areconstrainedtobenon-negative. f(x,y)inthex-directionaroundpoint(0,y)T.
Equations(6)and(7)definethed-dimensionaltrans-
formpairforgoingfromdirectx-spacerepresentation
totransformI-spacerepresentationandvice-versa.Nu- 2.2. DealingwithRotationin2D
mericalcomputationalgorithmscantakegoodadvan-
tageofthedimensionalfactorisationthatisinherentin In two dimensions, the Hermite transform is partic-
these equations. To avoid burdening the notation, we ularly suited for dealing with rotations (Koenderink
limitourselvesheretothe2Dcasebutourobservations and Van Doorn, 1992; Martens, 1997). This is be-
are readily generalized to any number of dimensions cause the Hermite coefficients fi,j can be expressed
d. For the 2D case, the multi-indices I are written as as multiple partial derivatives of an invariant scalar
pairsofnon-negativeintegers I = (i, j).Achangeof function(seeEq.(29)inSection3).Tensortheoryal-
variableisusedsoastobringthecentreoftheGaus- lows us to deduce that the set of fi,j coefficient with
sianwindowtotheorigin√ξ = (0,0)T,thecoordinate rank p ≡ i + j shouldbehavelikesymmetricrank p
unitsareselectedsothatσ 2=1andtheindependent covarianttensors(Spain,1960)underisometrictrans-
spacepositionalvariableiswrittenas(x,y)T.Equation formations (pure rotations and/or mirror symmetries
(6)cannowbereformulatedas aroundthewindowcentre).Anotherfamiliarinstance
of rank p tensors is that of degree p monomials xiyj
(cid:1)
1 (withp =i+j).Suchmonomialsbehavelikesymmet-
fi,j = π e−x2−y2Hi(x)Hj(y) ric contravariant tensors under linear transformations
x∈R2
× f(x,y)dxdy of the coordinates. Both covariant and contravariant
(cid:1) tensorsleadtoidenticallineartransformationlawsfor
+∞
1
= √ e−y2 H (y)f (y)dy thepurerotationtransformationsthataredealtwithin
π j i
y=−∞ thissection(Spain,1960).Wenowdeterminethema-
trixA(p)whichallowstheconversionofthemonomials
where xiyj underrotation,whichisthesameasthatneeded
(cid:1) +∞ to convert the Hermite coefficients fi,j . Since there
f (y)= √1 e−x2H(x)f(x,y)dx. (8) are p+1elementsofrank p,A(p)isa(p+1)2square
i π i
x=−∞ matrix.Asintheprevioussection,webringthecentre
of the Gaussian window to the origin ξ = (0,0)T. If
In other words, one may first take the 1D Hermite the coordinate axes are rotated by an angle θ so that
transform fi(y) along the x-axis for each of the y- any(x,y)T vectorischangedto(x(cid:11),y(cid:11))T accordingto
coordinates and then take the 1D Hermite transform
thelaw:
alongthey-axisforeachvalueofitogetthesetofco-
(cid:5) (cid:6) (cid:5) (cid:6)(cid:5) (cid:6)
efficients fi,j. A symmetric procedure allows getting x(cid:11) cosθ sinθ x
theinversetransforminasimilarmanner.Toestimate y(cid:11) = −sinθ cosθ y , (9)
the algorithmic complexity of those procedures, as-
sumethatnumericalintegrationalongthex-ory-axes
makeuseofLdiscretizationpointsandthattheiandj then,inthenewreferential,anymonomial(x(cid:11))i(y(cid:11))p−i
indicestakevaluesintheinterval[0,M−1].Thenum- (with0<i < p)canbeexpressedas
berofoperationsneededforeachofthetransformand
itsinverseisoftheorderofL2M+LM2.So,ifNstands
(x(cid:11))i(y(cid:11))p−i
forthelargestofLandM,thecomputationloadisof
theorderof2N3andforanynumberofdimensionsdit =(x.cosθ +y.sinθ)i(−x.sinθ +y.cosθ)p−i
canbeeasilyseentobeoftheorderofd×Nd+1.Ifthe (cid:3)j=p
orthogonalbasiswasnotdimensionallyseparable,the ≡ A(p)xjyp−j. (10)
i,j
correspondingcomputationalloadwouldhavebeenof j=0
Scale-SpaceImageAnalysisbasedonHermitePolynomialsTheory 129
The coefficients of matrix A(p) are evaluated by ex- If only rotations have to be dealt with, the polar
pandingthemiddleexpressionofEq.(10)intermsof Gauss-Laguerre transform (Jacovitti and Neri, 2000;
the monomials (x)j(y)p−j. As already stated above, KoenderinkandVanDoorn,1992;Martens,1997)may
wecandeducehowHermitecoefficientsareconverted bemoreefficientthantheCartesianHermitetransform.
underrotationnamely Thetwotransformtypesarestronglyrelatedsothatgo-
ingfromonetotheotherisasimplematter(Martens,
(cid:3)j=p 1997). However, it is always recommendable to start
fi(cid:11),p−i = Ai(,pj)fj,p−j. (11) fromthecomputationallyefficient(dimensionallysep-
j=0 arable)Cartesianformandthendeducethepolarform
whenneeded.
Thematrixcoefficientstakeaparticularlysimpleform
wheni = 0andwheni = p (topandbottomlinesof
2.3. Mehler’sFormulaanditsApplications
A(p))sincetheycanthenbeexpressed(usingEq.(10))
toScale-Space
fromthebinomialtheorem
(cid:5) (cid:6)
TodealwithHermitedomaintranslation,scalechange
A(p) =(−1)j p sinjθ.cosp−jθ and blurring, we suggest using a formula due to F.G.
0,j j
Mehler(late19thCentury,seeSzego¨ (1967),p.380).
Thisformulaisveryvaluablewhenstudyingthecon-
and
vergenceofHermitepolynomialexpansions.Mehler’s
(cid:5) (cid:6)
A(pp,)j = pj cosjθ.sinp−jθ (12) formula(r1e−adst2)−1/2exp(cid:5)2txy−t2(x2+y2)(cid:6)
(1−t2)
where (cid:3)∞ ti
(cid:5) (cid:6) (cid:5) (cid:6) = 2ii!Hi(x)Hi(y) (15)
p = p! = p (13) i=0
j (p− j)!j! p− j
foranyvalueoftintherange−1<t <1.Theproof
arebinomialcoefficients. ofthiskeyformulaisstraightforwardusing2DFourier
Moregenerally,thecoefficientsofmatrices A(p)can transform (Watson, 1933). B√y mult√iplying both sides
berecursivelycalculatedstartingwithtotalorder p = of Eq. (15) by e−x2 f(ξ +σ 2x)/ π where f is an
0. Since the zero order term is rotation invariant, the arbitrary 1D function and integrating both sides from
corresponding 1×1matrixhasitssingleentryequal −∞to+∞,onegets
tounityi.e. A(0) =1.UsingEqs.(9)and(10),weget (cid:1)
0,0 +∞ √
thefollowingrecursiverelations (cid:11) 1 e− (t1y−−tx2)2 f(ξ +σ 2x)dx
π(1−t2) x=−∞
Ai(,pj+1)=cosθ Ai(,pj)−sinθ Ai(,pj)−1 fori ≤ p (cid:3)∞ f (σ ; ξ)
= i tiH(y) (16)
2ii! i
and i=0
(cid:12)
A(pp++11,)j=cosθA(pp,)j−1+sinθA(pp,)j (14) where fi(σ; ξ) = √21πσ2 −+∞∞e−(x2−σξ2)2Hi(σx−√ξ2)
f(x) dx is the Hermite coefficient of f for a Gaus-
sian window centred at x = ξ and of size σ. With
with the convention that the matrix entries A(p) with changeofvariables,thisresultcanberewrittenas
i,j (cid:11)
j > 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 prod-
uct 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 pro-
cessing, 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).