domain: math spi phys

text:
Wenxiang Ding
email: wenxiang.ding@insa-cvl.fr

Maxime Bavencoffe
email: maxime.bavencoffe@insa-cvl.fr

Marc Lethiecq
email: marc.lethiecq@insa-cvl.fr

Accurate Coupled Vibration Analysis of a Piezoelectric Ceramic Cylinder by the Superposition Method

Keywords: coupled vibration, analytical model, superposition method, electrical impedance, piezoelectric cylinder, finite element method

The radial and thickness extensional vibration modes in piezoelectric cylinders are always inevitably coupled due to the finite dimension, Poisson's ratio, and piezoelectric effect. In this paper, an analytical model based on the superposition method is developed to obtain the coupled dynamic response of a piezoelectric cylinder under an applied voltage. The problem can be described mathematically by three partial differential equations with mixed boundary conditions in the cylindrical coordinates system. To solve this, the problem is decomposed first into two building blockvibrations in radial and thickness directions. In each building block, the expressions of displacements and electric potential are assumed and then the induced dynamic responses, such as in-plane stress and electric displacements, are calculated. Finally, the vibration responses of the two building blocks are superimposed to satisfy the mixed boundary conditions using Fourier and Fourier-Bessel series expansions. Electrical impedance of a typical piezoelectric disk and frequency spectrum of piezoelectric cylinders of different diameter-to-thickness ratios are calculated by the present analytical method as well as by finite element method. Comparison shows an excellent agreement. This analytical model can be applied to material characterization and the design and the optimization of the active elements of piezoelectric devices.

Introduction

Piezoceramic cylinders of different diameter-to-thickness (D/T) ratio have been widely used as the active element of ultrasonic transducers in many applications such as medical diagnostics, therapy, non-destructive evaluation (NDE), cleaning, atomizing liquids, drilling or milling materials, welding plastics as well as for many other purposes [START_REF] Farmer | The application of power ultrasound to the surface cleaning of silica and heavy mineral sands[END_REF][START_REF] Marechal | High-frequency transducers based on integrated piezoelectric thick films for medical imaging[END_REF][START_REF] Lethiecq | Piezoelectric transducer design for medical diagnosis and NDE[END_REF][START_REF] Bavencoffe | Modelling and experimental measurements of the mechanical response of piezoelectric structures from millimetre to micrometre[END_REF]. Depending on the material properties, the dimensions, the electric field applied and the poling direction, two principal vibration modes can be excited: the thickness extensional and the radial modes. One dimensional (1D) analytical models such as Van Dyke [START_REF] Dyke | The Piezo-Electric Resonator and Its Equivalent Network[END_REF], Mason [START_REF] Mason | Electromechanical transducers and wave filters[END_REF], or KLM [START_REF] Krimholtz | New equivalent circuits for elementary piezoelectric transducers[END_REF] models have been proposed to characterize the dynamic response of the resonator when the diameter is much larger (i.e. ceramic disks) or much smaller (i.e. ceramic rods) than the thickness. However, these 1D models are not applicable to describe the coupling between radial and thickness modes. To address this, a two-dimensional (2D) or three-dimensional (3D) model is needed.

Due to its axisymmetric structure, the 3D vibration analysis of a piezoceramic cylinder can be reduced to a 2D problem in the cylindrical coordinate and can be described mathematically by a system of three partial differential equations (PDEs). The biggest difficulty of solving this problem lies in finding the appropriate forms for the unknown variables, i.e. mechanical displacements and electric potential, to satisfy both the differential equations and boundary conditions. Different kinds of potential functions were proposed in literature by Saito [START_REF] Saito | The Axially Symmetrical Deformation of a Short Circular Cylinder[END_REF], Rajapakse et al [START_REF] Rajapakse | Electroelastic field of a piezoelectric annular finite cylinder[END_REF][START_REF] Chen | Electric Charge Loading of a Piezoelectric Solid Cylinder[END_REF], Wei and Chau [START_REF] Wei | Analytic Solution for Axial Point Load Strength Test on Solid Circular Cylinders[END_REF][START_REF] Wei | Analytic Solution for Finite Transversely Isotropic Circular Cylinders under the Axial Point Load Test[END_REF], and Wang et al [START_REF] Wang | A general solution and the application of space axisymmetric problem in piezoelectric material[END_REF]. In these works, one or two potential functions were introduced, then the mechanical displacements and electric potential could be represented by a differential form to satisfy the equations of equilibrium and boundary conditions. However, they are only suitable in static regimes, not in dynamics. Recently, Ding et al [START_REF] Ding | General solutions for coupled equations for piezoelectric media[END_REF] obtained the general solutions for coupled dynamic equations for a transversely isotropic piezoelectric medium by introducing two potential functions. Although the method was later extended to calculate the fundamental frequencies of a circular plate under open and closed circuit conditions [START_REF] Ding | Exact solutions for free vibration of transversely isotropic piezoelectric circular plates[END_REF], it is limited to the two specific boundary conditions, named as elastic simple support and rigid slipping support.

Approximate methods by virtue of weakening the coupling effect were developed lately by Brissaud [START_REF] Brissaud | Characterization of piezoceramics[END_REF][START_REF] Brissaud | Three-dimensional modeling of piezoelectric materials[END_REF], Iula [START_REF] Iula | An approximated 3-D model of cylindershaped piezoceramic elements for transducer design[END_REF], and Zhang [START_REF] Zhang | Three-dimensional coupled vibration theory for the longitudinally polarized piezoelectric ceramic tube[END_REF]. The displacement along each propagation direction was considered as independent, i.e. only related to the corresponding coordinate, and the electrical quantities in non-polarized directions were assumed to be zero. The shear deformation was always neglected. In some ways, these assumptions did not satisfy the constitutive equations. A set of 2D second order approximate equations were deduced by Lee [START_REF] Lee | A new two-dimensional theory for vibrations of piezoelectric crystal plates with electroded faces[END_REF][START_REF] Huang | Extensional, thickness-stretch and symmetric thickness-shear vibrations of piezoceramic disks[END_REF][START_REF] Lee | An accurate two-dimensional theory of vibrations of isotropic, elastic plates[END_REF]. The calculated dispersion curves and resonance frequencies of piezoelectric crystal plates gave a good prediction. Similar second order equations were also developed by Li [START_REF] Li | Two-dimensional equations for piezoelectric thin-film acoustic wave resonators[END_REF][START_REF] Li | Forced coupling vibration analysis of FBAR based on two-dimensional equations associated with state-vector approach[END_REF][START_REF] Li | Two-Dimensional Plate Theory for the Analysis of Coupling Vibrations in Shear Mode FBARs[END_REF] and were applied to the vibration analysis of thin-film bulk acoustic wave resonators.

The coupling between vibration modes has been widely studied. Satisfactory results have been obtained in static regime for rectangular or cylindrical geometries. To the best of the author's knowledge, the dynamic coupled response of a piezoelectric ceramic under external electric field has not been established. In the present work, a complete analytical solution for coupled vibration in piezoelectric ceramic cylinders in vacuum is proposed based on the superposition method, which was first proposed by Gorman [START_REF] Gorman | Free vibration analysis of the completely free rectangular plate by the method of superposition[END_REF][START_REF] Gorman | Free in-plane vibration analysis of rectangular plates by the method of superposition[END_REF][START_REF] Gorman | Accurate in-plane free vibration analysis of rectangular orthotropic plates[END_REF]. In section 2, the mathematical formulation is presented. Section 3 shows the coupling effect between radial and thickness modes and compares the results with those calculated by finite element (FE) method. Conclusions are given in section 4.

Mathematical Formulation

Here, a piezoelectric cylinder of thickness 2H and diameter 2R in vacuum is considered. Fig. 1 shows the coordinate system and dimensions. The poling direction is parallel to z-axis. The cartesian coordinates are x, y, and z. They correspond to r, θ, and z coordinates in a cylindrical system. An AC-voltage is applied between the top and bottom surfaces covered by metallic electrodes. Due to the axisymmetric structure, the 3D vibration analysis can be reduced to a 2D axisymmetric one as shown in Fig. 1 The constitutive equations of a piezoelectric element can be expressed as [START_REF] Royer | Elastic waves in solids I: Free and guided propagation[END_REF] , 1, 2,3 and , 1,..., 6

E i i S i ij j i T c S e E i j D E e S α αβ β α α α α β ε = - = = = + (1)
where E c αβ are the elastic stiffness constants under constant electric field, i e α are the piezoelectric constants,

S ij

ε are the dielectric constants under constant strain, T α , S β , i E and i D are the stress, strain, electric field and electric displacement, respectively.

For the 2D axisymmetric model, the interaction between θ and r/z-axis is omitted, thus Srθ = Sθz = 0 and Trθ = Tθz = 0. Based on this, the constitutive equations can be simplified and rewritten in the cylindrical coordinate system as follows

= = = = D = ε D = ε E E E r r z z E E E r z z E E E z r z z E rz rz r S r
For axisymmetric structures, the strain and displacements are related by , , ,

r z rz u u w u w S S S S r r z z r θ ∂ ∂ ∂ ∂ = = = = + ∂ ∂ ∂ ∂ (3)
The electric fields are related to electric potential by ,

r z E E r z ϕ ϕ ∂ ∂ = - = - ∂ ∂ (4)
The coupled motion equations and Maxwell's equation for electrostatic field are

T T T T u r z r t T T T w r z r t D D D r z r θ ρ ρ - ∂ ∂ ∂ + + = ∂ ∂ ∂ ∂ ∂ ∂ + + = ∂ ∂ ∂ ∂ ∂ + + = ∂ ∂ ( 5
)
where ρ is the mass density of the material. Substituting Eqs. ( 2)-(4) into Eq. ( 5) yields three coupled differential equations

( ) ( ) ( ) ( )
ϕ ρ ϕ ϕ ρ ϕ ϕ ε ε ∂ ∂ ∂ ∂ ∇ - + + + + + = ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ + ∆ + ∆ - + + ∆ + = ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ + ∆ + ∆ + -∆ - = ∂ ∂ ∂ ∂ ∂ (6)
where u, w are the displacements along r and z axes and φ is the electric potential. The operators ∆ and ∇ 2 are defined as

2 2 2 2 1 1 1 , r r r r r r r ∂ ∂ ∂ ∂ ∆ = + ∇ = ∆ = + - ∂ ∂ ∂ ∂ (7)
Since there are no electrodes on the lateral surface, the electric displacement is equal to zero. In addition, the piezoelectric element is stress free and connected to an ACvoltage of angular frequency ω, therefore, the boundary conditions can be expressed as 0 0 at 0, at

r rz r j t z rz T T D r R T T e z H ω ϕ ϕ = = = = = = = ± = ± (8)
Owing to the piezoelectric effect, a symmetric constraint is imposed on the sample when an AC-voltage is applied. This leads to the displacement u(r,z) being an even function versus z coordinate while w(r,z) and φ(r,z) are odd functions of z. It can be summarized as follows

( ) ( ) ( ) ( ) ( ) ( ) , , , , , , u r z u r z w r z w r z r z r z ϕ ϕ = - = - - = - - (9)

General solution for the vibration responses

Since the piezoelectric element is axisymmetric with respect to z-axis and symmetric with respect to r-axis, a half part is chosen to be studied, as shown in Fig. 2(a). Two pairs of small circles indicate the symmetry condition along each side. As mentioned above, the major difficulty of solving the PDEs is to find an appropriate expression of the variables such that the PDEs and boundary conditions can be satisfied simultaneously. To address this, the model is divided into two building blocks as shown schematically in Fig. 2(b). In fact, the superposition method adopted here is a modified version of the method of separation of variables in mathematics. This separation process is tightly related to normal vibration modes in radial and thickness directions. Let's now consider the first building block on the right of Fig. 2(b). Supposing that only the radial vibration modes are induced, the displacements and electric potential are expressed as follows [START_REF] Chen | Electric Charge Loading of a Piezoelectric Solid Cylinder[END_REF][START_REF] Senjuntichai | Piezoelectric cylinder under voltage and axial loading[END_REF] ( )

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 0 0 , , ( 1, 2,
..., + ) , m m m m m m m m m u r z U z J k r w r z W z J k r m r z z J k r ϕ ψ = = = ∞ = ( 10
)
where J0(kmr) and J1(kmr) are Bessel functions of the first kind of order zero and one, respectively. The wavenumber in the radial direction, km, is the m-th root of J0(kmR) = 0. For simplicity, the time component e jωt is omitted.

By substituting Eq. ( 10) into the coupled PDEs Eq. ( 6), we obtain ρω =0

( ) ( ) ( ) ( ) ( ) ( ) ( )
E E E E m m m m m m m c U z c k U z c c k W z k e e z ψ ′′ ′ ′ + - - + - + (11a) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2
( ) ( ) ( ) ( ) 1 2 3 4 0 m m m m m m m m a U z a U z a W z a W z ′′ ′′′ ′ + + + = (12a) ( ) ( ) ( ) ( ) 1 2 3 4 0 m m m m m m m m b U z b U z b W z b W z ′′′ ′ ′′ + + + = (12b)
where the coefficients ami and bmi (i = 1,2,3,4) are constants related to material properties, ω, and km. They are listed in Appendix (Eq. (A. 1)). Differentiating Eq. (12b) once with respect to z and combining with Eq. (12a), the quantities

( ) m W z ′′′ and ( ) m W z ′ can be expressed in terms of ( ) iv m U z , ( ) m U z ′′ , and
( ) m U z . Let ( ) m W z ′′′ be equal to ( ) m W z ′
differentiating twice with respect to z, a sixth order homogeneous differential equation involving only

( ) m U z and its derivatives is obtained ( ) ( ) ( ) ( ) 1 1 1 0 vi iv m m m m U z bU z cU z dU z ′′ + + + = (13)
where ( ) ( )

1 3 2 4 1 1 3 3 1 1 4 2 1 4 2 3 3 1 1 2 4 3 1 , m m m m m m m m m m m m m m m m m m m m b a b a b a b a b c a b a b a b a b d a b a b = + - = - - = -
Substituting the exponential function z e λ into Eq. ( 13), the characteristic equation related to the quantity

( ) m U z is obtained 6 4 2 1 1 1 0 b c d λ λ λ + + + = (14)
For given ω and km, three roots with respect to λ 2 are given by the roots of a general cubic equation [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF][START_REF] Weisstein | Cubic formula[END_REF]. The roots of the characteristic equation are represented by three sets of square roots, denoted as ( )

1, 2,3 i mi k i λ = ± = (15)
The frequency ω and wavenumber km are real, but

2 i
λ may be either real-positive, real-negative, or complex, which leads to several cases for the roots kmi. In the present work, roots are considered as complex for the sake of generalization. Although we may have + kmi and -kmi results, only the + kmi need to be considered since they give the same solution after the conversion from exponential functions to hyperbolic ones [START_REF] Tiersten | Linear Piezoelectric Plate Vibrations[END_REF].

The general solution ( ) m U z , being an even function of z, can be expressed as

( ) ( ) 3 1 cosh m mi mi i U z U k z = = ∑ (16)
where Umi are arbitrary constants.

Turning back to Eqs. ( 11) and ( 12), the quantities

( ) m W z and ( ) m z ψ can be expressed by ( ) m U z ( ) ( ) ( ) ( ) ( ) 3 2 2 2 1 sinh v m m m m mi mi mi i W z bU z c U z d U z B U k z = ′′′ ′ = + + = ∑ (17) ( ) ( ) ( ) ( ) ( ) 3 3 3 3 1 sinh v m m m m mi mi mi i z bU z cU z d U z C U k z ψ = ′′′ ′ = + + = ∑ (18) where ( ) ( ) ( ) ( ) ( ) ( )
= + + = - - + = - + - - = - - - = + + = - - = - ( )( ) ( ) ( ) (
)

1
+ - - - = - + - -
and the coefficients cmi (i = 1,2,3) are given in Appendix (Eq. (A. 2)). Now let's turn to the second building block on the right of Fig. 2(b). The derivation is similar to the one of the first building block but slightly more complicated because of the cylindrical coordinate system. Therefore, for completeness sake, the derivation process is briefly presented here. In the second building block, the thickness vibrations are induced. To avoid confusion, a n-index is used for the notations of the displacements and electric potential [START_REF] Gorman | Free in-plane vibration analysis of rectangular plates by the method of superposition[END_REF][START_REF] Kim | Coupled vibration analysis for a piezoelectric array element using superposition method[END_REF]

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) , cos , sin and 1, 2,..., , n 2 1 si 2 n n n n n n n n n n u U r r k n r z k z w r z W k k H r z n r z z π ϕ ψ = = = +∞ = - = (19)
Inserting this into Eq. ( 6), we have ρω 0

( ) ( ) ( ) ( ) ( ) ( ) ( ) 2
E E E E n n n n n n n c U r c k U r c c k W r e e k r ψ ′ ′ ∇ + - + + + + = (20a) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 2 13
55 55 33 15 33 ρω 0

E E E E n n n n n n n n c c k U r c W r c k W r e r e k r ψ ψ ′ ′ - + ∆ + ∆ + - + ∆ - = (20b) ( ) ( ) ( ) ( ) ( ) ( ) 2
′ -+ ∆ + ∆ - -∆ + = (20c)
where the apostrophe implies differentiation with respect to variable r. By combining the equations following a procedure similar to the previous case, we can eliminate the quantity ( )

n r ψ . Then we have ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 2 1 2 3 4 2 1 2 3 4 0 0 n n n n n n n n n n n n n n n n a U r a U r a W r a W r b U r b U r b W r b W r ′ ′ ∇ + + ∇ + = ′ ∆∇ + ∆ + ∆ + = (21)
After combining these equations, we can eliminate the quantity

( ) ( ) ( ) ( ) 6 4 2 1 1 1 0 n n n n U r b U r c U r dU r ∇ + ∇ + ∇ + = (22)
where the definition of the operator 2 ∇ is given in Eq. ( 7).

For given ω and kn, the three roots for ( )

2 n U r ∇ are denoted as 2 ni
k , then, the solution of Eq. ( 22) is given by ( ) ( )

3 1 1 n ni ni i U r U I k r = = ∑ (23)
where Uni are arbitrary constants and ( )

1
( ) ( ) ( ) ( ) ( ) 3 4 2 2 2 2 0 1 n n n n ni ni ni i W r b U r c U r d U r B U I k r = = ∆∇ + ∆∇ + ∆ = ∑ (24) ( ) ( ) ( ) ( ) ( ) 3 4 2 3 3 3 0 1 n n n n ni ni ni i r b U r c U r d U r C U I k r ψ = = ∆∇ + ∆∇ + ∆ = ∑ (25)
where ( )

0 ni
I k r is the modified Bessel function of the first kind of order zero.

The coefficients ani, bni and cni (i = 1,2,3,4) for the second building block are detailed in the Appendix (Eq. (A. 3)). The other coefficients, such as bi, ci and di (i = 1,2,3), are exactly the same as those that appeared in the first building block. To obtain these coefficients, we need only to replace symbol m by n. For clarity, they are also given in the Appendix (Eq. (A. 4)).

With the solutions of the two building blocks above, the general solution for mechanical displacements and electric potential can be expressed as

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 1 3 1 1 1 1 1 , + c = os cos c h os m m m n mi mi m ni n n n m i n n i u r z U z J k r r U k z J k r U I k r U k z k z ∞ ∞ = = ∞ = = ∞ =      = +  ∑ ∑ ∑ ∑ ∑ ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 0 1 1 1 1 0 0 1 3 s , + sinh sin in m m m n mi mi mi m n n n n n i m i ni i n w r W k z k r z W z J k r r B U k k z z J r B U I k ∞ = ∞ = = ∞ ∞ = =  =  =   +   ∑ ∑ ∑ ∑ ∑ (26) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 0 3 1 0 0 1 1 1 1 , + sinh si i s n n n n m m m m n mi mi mi ni ni ni n n i m r z z J k r r C U k z J k r k k z k C U I r z ϕ ψ ψ ∞ ∞ ∞ = = = = ∞ = = +   =     ∑ ∑ ∑ ∑ ∑
The unknown coefficients mi U and ni U can be obtained after substituting the solution Eq. ( 26) into the boundary conditions Eq. [START_REF] Saito | The Axially Symmetrical Deformation of a Short Circular Cylinder[END_REF]. Once these coefficients are determined, the following expressions can be used to evaluate the dynamic vibration responses of stress, electric field, and electric displacements inside the element: , cosh cos , , sinh sin

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 3
m E E E E r mi mi mi m m mi mi i m ni E E E E ni ni n ni ni n ni i n E E mi mi mi m m J k r T r z c B e C k c k J k r c c k z U r I k r c B e C k c k I k r c c k z U r T r z c B e C k c k J k r c θ ∞ = = ∞ = =     = + + + -           + + + + -         = + + +   ∑∑ ∑∑ ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 3
m E E mi mi i m ni E E E E ni ni n ni ni n ni i n E E z mi mi mi m m mi mi i m E E ni ni n J k r c k z U r I k r c B e C k c k I k r c c k z U r T r z c B e C k c k J k r k z U c B e C k c ∞ = = ∞ = = ∞ = =   -         + + + + -         = + +   + + + ∑∑ ∑∑ ∑∑ ( ) ( ) 3 0 1 1 cos ni ni n ni i n k I k r k z U ∞ = =     ∑∑ ( ) ( ) ( ) ( ) ( ) ( ) ( ) 3
E E rz mi mi m mi m mi mi i m E E ni ni ni n ni n ni i n T r z c B e C k c k J k r k z U c B e C k c k I k r k z U ∞ = = ∞ = =   = - - +     + + -   ∑∑ ∑∑ (27) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 3 3 1 1 1 1 1 1 3 3 0 0 1 1 1 1 , sinh sin , cosh cos r mi m mi m mi ni ni ni n ni i m i n z mi mi mi m mi ni n ni n ni i m i n E r z C k k z J k r U C k I k r k z U E r z C k k z J k r U C k I k r k z U ∞ ∞ = = = = ∞ ∞ = = = = = - = - - ∑∑ ∑∑ ∑∑ ∑∑ ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 3
ε sin D , ε cosh ε S r mi mi m mi mi m mi i m S ni ni ni n n ni ni i n S z mi mi mi m mi m mi i m S ni ni n ni r z e B C k e k k z J k r U e B C k e k k z I k r U r z e B C k e k k z J k r U e B C k e k ∞ = = ∞ = = ∞ = =   = - +     + - -     = - +     + - +   ∑∑ ∑∑ ∑∑ ( ) ( ) 3 0 1 1 cos n ni ni i n k z I k r U ∞ = =

∑∑

Contrary to the common assumption [START_REF] Brissaud | Characterization of piezoceramics[END_REF][START_REF] Iula | An approximated 3-D model of cylindershaped piezoceramic elements for transducer design[END_REF][START_REF] Zhang | Three-dimensional coupled vibration theory for the longitudinally polarized piezoelectric ceramic tube[END_REF], the electric filed in radial direction Er is not equal to zero. In addition, the symmetry property versus z-axis, shown in Eq. ( 9), can be verified here.

Determination of unknown coefficients

Assuming the number of vibration modes induced in the radial and thickness directions are K1 (m = 1, 2,…,K1) and K2 (n = 1,2,…,K2), there will be 3(K1+K2) linear equations for the 3(K1+K2) unknown coefficients of Umi and Uni by substituting the solutions into the boundary conditions of Eq. ( 8). The 3(K1+K2) linear equations can be presented in a concise way as

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 2 1 1 2 2 1 3 1 0 1 1 3 1 1 1 1 1 3 2 2 1 1 1 1 3 2 2 1 1 1 0 as 0 cos 0 as 0 0 as 0 sin 0 as K mi mi m z m i K K ni mi ni mi n r n i m K K mi ni mi ni m rz m i n K K ni mi ni mi n rz n i m F U J k r T z H F U G U k z T r R F U G U J k r T z H F U G U k z T r R = = = = = = = = = = = = = ± =   + = = =       + = = ± =       + = =     ∑∑ ∑∑ ∑ ∑∑ ∑ ∑∑ ∑ ( ) ( ) ( ) ( ) 1 2 2 1 3 3 3 0 0 0 1 1 1 3 3 3 1 1 1 0 as sin 0 as 0 K K mi ni mi mi ni m m i n K K ni mi ni mi n r n i m F U G U J k r z H F U G U k z D r R ϕ ϕ ϕ = = = = = = =   + = = ± = ±       + = = =     ∑∑ ∑ ∑∑ ∑ (28)
where mi j F and ni j F (j =1,2,3) can be easily obtained from Eq.( 27); they are given by ( )

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1
F c B C e k c k k H F c B C e k c k I k R c c I k R R F c B e C k c k k H F c B e C k c k I k R F C k H F   = + +     = + + +     = - - +     = + -   = = - ( ) ( ) 11 15 15 1 ε 1, S ni ni ni n ni C e B k e k I k R   + -   (29)
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1
H mi E E m n mi R ni E E ni ni ni n n ni m m H mi E E mi mi m mi m mi n R ni ni n ni m m G c c J k R k z k z dz RH G c B e C k c k k H rI k r J k r dr J k r G c B e C k c k J k R k z k z dz H G C k H rI k r J k r dr J k =   = + -     = - - +   = ∫ ∫ ∫ ∫ ( ) ( ) ( ) ( ) ( ) 2 3 11 15 15 1 0 ε 2 sinh sin H mi S mi mi m mi m mi n r G C e B k e k J k R k z k z dz H   = - +   ∫ (30)
Eq. [START_REF] Gorman | Accurate in-plane free vibration analysis of rectangular orthotropic plates[END_REF] shows that the stress z T , rz T , and electric potential ϕ at z H = ± can be expanded with the first kinds of Bessel functions of order 0, 1, and 0, respectively. The stress r T , rz T , and electric displacement r D at r R = can be expanded with cosine, sine, and sine functions, respectively. The boundary conditions are satisfied in the form of Fourier series and Fourier-Bessel series expansions. The summation of the expanded coefficients should be equal to zero except for the electric potential at z H = ± . The nonhomogeneous linear system in matrix notation can be described by = Ax b [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF] where A is a 3(K1+K2)×3(K1+K2) coefficient matrix, x is the unknown variable of Umi and Uni, and b is a 3(K1+K2)×1 vector with K1 non-zero elements related to φ0. For clarity, an example is shown here, where K1 = 2 and K2 = 3. The augmented matrix of the linear system is assembled and depicted schematically in Fig. 3

Electrical impedance

The electrical impedance of piezoelectric material is defined as

= V Z I ( 32
)
where V is the voltage applied between the electrodes and I is the electric current flowing through the sample.

Since the electric potential on the top and bottom electrodes are

0 j t e ω ϕ ± , we have 0 2 H H j t z H H V E dz dz e z ω ϕ ϕ - - ∂ = = - = - ∂ ∫ ∫ (33)
The electric current is defined by the differential of charge versus time. The charge on electrodes can be calculated by integrating the electric displacement density D3 over the surface. Therefore,

( ) ( ) 3 0
, , 2 ,

R j t z S I D r H t dS j e D r H rdr t

ω π ω ∂ = = ∂ ∫ ∫ (34)
Due to Eq. ( 19) and Eq. ( 27), ( )

, z D r H is given as ( ) ( ) ( ) ( ) 3 33 33 31 0 1 1 D , ε cosh S z mi mi mi m mi m mi i m r H e B C k e k k H J k r U ∞ = =   = - +   ∑∑ ( 35
)
After obtaining all the unknown coefficients in Section 2.2, we can substitute them into Eqs. ( 32)-( 35), then the electrical impedance is determined. The time component j t e ω and the constant φ0 vanish in the results.

Results and Discussion

To verify the validity of the analytical model proposed here, the results are compared with those computed by a FE method. The FE analysis is carried out with COMSOL Multiphysics software in vacuum. Quadratic Lagrange rectangular elements with nine nodes are used and a frequency domain analysis is performed using 8 elements per shear wavelength for the highest frequency. The validation is performed in two parts: the electrical impedance and mode shape, represented by the spatial distribution of the mechanical displacements and electric potential, given in section 3.1. The predicted frequency spectrum comparison of different diameter-to-thickness ratios is shown in section 3.2.

Soft PZT materials provide high coupling factors and piezoelectric charge coefficients. They are widely used in applications such as medical diagnostic transducers, doppler flow meters, accelerometers and NDE. A soft PZ27 (Ferroperm Piezoceramics [START_REF]Ferroperm piezoceramic material data for modelling[END_REF]) ceramic disk is used here, with diameter of 16 mm and thickness of 1.13 mm. An AC-voltage of ±1V is applied between the top and bottom electrodes.

Material properties are listed in

Electrical impedance and mode shape of a ceramic disk

Since the analytical model is based on Fourier and Fourier-Bessel series expansions, any desired degree of accuracy could be achieved by increasing the number of terms K1 and K2. To determine the values of K1 and K2, a convergence test for resonance frequency (fr) is performed. Fig. 4 shows the variation of resonance frequency with number of terms K1, for all vibration modes in the frequency range between 0 and 2.4 MHz. The number of terms K2 is set at 2 and the frequency resolution is set at 1kHz.

One can see that as the number of terms K1 increases, more and more radial or coupling modes appear -in other words, to characterize high order radial harmonics and the coupling between thickness and radial modes in higher frequency ranges, K1 needs to increase. One can also see that the calculated values of fr converge vertically more or less to an asymptotic value as K1 increases. The upper left area in Fig. 4, separated by a symbolic oblique solid line, can be considered as the area of convergence. It should be noted that in the case presented here (from 0 MHz to 2.4 MHz), 21 K1 terms and 2 K2 terms are sufficient to allow this asymptotic value to be obtained with good reproducibility. 5. To measure the discrepancy, a mean absolute percentage deviation (MAPD) metric [START_REF] Wang | Monitoring of the strength gain of concrete using embedded PZT impedance transducer[END_REF] is introduced and defined as

( ) ( ) ( ) 1 1 MAPD= 100 N Analytical i FE i i FE i Z f Z f N Z f = - × ∑ (36)
where Z(fi) is the electrical impedance value at frequency fi and N is the number of points evaluated. The agreement between the analytical and the FE model improves and the MAPD decreases from Fig. 5(a) to Fig. 5(c). The MAPDs of the three sets of results are close to 55%, 15% and 2%, respectively. For a better illustration, the resonance and antiresonance frequencies in the third case are extracted and compared with those of the FE model. They are listed in Table 2, including the discrepancies in parentheses. The values are very close: the maximum discrepancy is inferior to 3 kHz for the frequency range calculated and the maximum deviation is under 2%, which occurs at the lower frequencies. , and electric potential φ are calculated by Eq. ( 26). Results at the antiresonance frequency of the fundamental radial (R1) and thickness (T1) modes, as labelled in Fig. 5, are shown in Fig. 6 and Fig. 7. Comparison shows an excellent agreement, in accordance with [START_REF] Kocbach | Finite element modeling of ultrasonic piezoelectric transducers[END_REF], including not only the spatial distribution but also the amplitudes (i.e. total displacement and electric potential respectively). The vibration pattern of R1 mode can be easily recognized from Fig. 6(a) & (c): the piezoceramic disk extends in the radial direction while keeping close to constant in the thickness direction. The electric potential distribution, as shown in Fig. 6(b) & (d), is uniform in the radial direction but linear in its poling direction, i.e. the thickness direction. From Fig. 7(a) & (c) we can see that the vibration of T1 mode is predominant in the thickness direction but is greatly influenced by the higher harmonics of radial modes. Additionally, the displacement is higher at the surface center than at the periphery. Owing to the piezoelectric effect, the distribution of electric potential of T1 mode (Fig. 7(b) &(d)) is similar to that of displacement. Besides, although an ACvoltage of ±1V is applied between the top and bottom electrodes in the poling direction, the electric potential inside some locations of the element exceeds this range and is more than ±2V.

In addition, compared to the FE method, the analytical method shows a great advantage in terms of computational efficiency: the calculation time is reduced from the FE method's 181.3s to only 2.8s for the third case. Thus, a comparable accuracy is attained 64 times faster.

Frequency spectrum

The present analytical model has been extended to predict the frequency spectrum of PZ27 piezoelectric disks, that is, relations between resonance frequencies and diameterto-thickness (D/T) ratios, as shown in Fig. 8. Here, the material is assumed to be lossless. The D/T ratio varies from 0.2 to 20 with a step of 0.05 while the thickness is kept constant as 1.13 mm. This range covers most of the piezoelectric element geometries used in ultrasonic transducers. Since the electrical impedance is inversely proportional to the electrode area, the amplitude is first normalized to a unit area and then normalized between 0 and 1. The resonance amplitude of the vibration modes is represented by the color scale in Fig. 8. Besides, to avoid a convergence analysis in each case, K1 and K2 are set at 30 and 2 for D/T ratios larger than 5 and are both set at 10 for the other cases.

The overall vibration pattern shown in Fig. 8 is similar to the measured or calculated frequency spectrum by FE model [START_REF] Kocbach | Finite element modeling of ultrasonic piezoelectric transducers[END_REF][START_REF] Ikegami | Frequency spectra of resonant vibration in disk plates of PbTiO3 piezoelectric ceramics[END_REF][START_REF] Guo | Measurement and prediction of the frequency spectrum of piezoelectric disks by modal analysis[END_REF][START_REF] Dimi | Resonance frequencies of PZT piezoceramic disks: a numerical approach[END_REF]. Results from the FE model, marked by black circles, are superimposed over analytical ones in Fig. 8. An excellent agreement can be observed. Several features are also clearly shown in the frequency spectrum. The three "terraces", as mentioned in [START_REF] Guo | Measurement and prediction of the frequency spectrum of piezoelectric disks by modal analysis[END_REF], are noted as edge modes (E), thickness shear modes (TS), and thickness extensional modes (TE). Their locations are also indicated by the arrows on the right of the plot area. The TS terrace is the most obvious in Fig. 8. It is close to flat and the number of TS modes increases as D/T ratio increases. The TE terrace is obviously characterized by its high resonance amplitude, i.e. the lighter color. The E terrace can be speculated because the rate of change of resonance frequency with D/T ratio is reduced in this region. The two horizontal dotted lines, labeled as k33 and kt in Fig. 8, are results from the one-dimensional (1D) KLM model [START_REF] Krimholtz | New equivalent circuits for elementary piezoelectric transducers[END_REF]. As D/T ratio approaches 0, the piezoelectric cylinder turns into a long rod with small section poled along its largest dimension, which can be modelled by KLM model with electromechanical coupling coefficient equal to k33 (D/T → 0). As D/T ratio approaches 20 or even higher, the cylinder turns into a thin disk poled along its thickness. It can be modelled by KLM with coupling coefficient kt (D/T → +∞). The electrical impedances of these two extreme cases calculated by two analytical models are shown in Fig. 9. In Fig. 9(a), only the length expander (LE) mode exists in the frequency range (0 to 2.4 MHz), and the resonance frequency is close to 1.1 MHz. In Fig. 9(b), the TE mode is highly perturbated by other modes, and the resonance frequency is around 1.7 MHz. These two resonance frequencies, calculated by model, could be considered as two asymptotic bounds for real situation: an upper bound (D/T → 0) and a lower bound (D/T → +∞). They are found by our analytical model with a discrepancy of less than 1%. This small discrepancy is linked to the D/T ratio: contrary to KLM model, the radial contribution is taken into account in our analytical model by means of a finite value of D/T. Plus, compared with KLM model, the analytical model proposed here can describe not only the two fundamental vibrations in radial and thickness directions but also the coupled modes such as E and TS modes.

Conclusions

In this paper, accurate analytical solutions have been obtained for coupled vibrations in piezoceramic cylinders based on superposition method. The problem has been divided into two building blocks -vibrations in radial and thickness directions -and solution in each building block are superimposed to form the final solution. The mixed boundary conditions are satisfied using Fourier and Fourier-Bessel series expansions.

The solutions for the in-plane dynamic responses, such as stress, electric field, and electric displacements, are formulated. Electrical impedance and mode shape of a typical size PZ27 piezoelectric disk, as well as the frequency spectrum of the same material but with different diameter-to-thickness ratios are calculated and compared to those of the FE method. An excellent agreement is observed. The mean absolute percentage deviation (MAPD) of the electrical impedance curve is inferior to 2%. In addition, the analytical method shows a great advantage in terms of computational efficiency -it is over 64 times faster than FE one for the same level of precision. Furthermore, comparison shows that the present analytical method is more accurate than KLM model: not only the radial and thickness modes but also the other coupled modes can be characterized using the analytical model proposed here.

The present analytical model can be used to help with the characterization of some material parameters but also for the design and the optimization of the active elements of piezoelectric devices, such as ultrasonic transducers.

( ) ( ) The constants cmi (i = 1,2,3) appearing in Eq. ( 18) are defined by ( The coefficients ani, bni and cni (i = 1,2,3,4) for the second building block (on the right of Fig. 2(b))) are defined as (

Fig. 1 .

 1 Schematic representation of the piezoelectric cylinder element. (a) threedimensional (3D) model; (b) two-dimensional (2D) axisymmetric model.

  The radial km and thickness kn wavenumbers of vibration modes are induced in the first and second building blocks respectively (Fig 2.(b)). After that, their vibration responses are superimposed to satisfy the boundary conditions. (a) (b) Fig. 2. Schematic representation of a half of the 2D axisymmetric piezoelectric element (a) and two divided building blocks (b) for the general solution.

  ni I k r is the modified Bessel function of the first kind of order one. The solutions for

  =1,2,3) are obtained using the Projection method[START_REF] Kevorkian | An accurate method for free vibration analysis of structures with application to plates[END_REF]

  the rest are zeros. Once the constant φ0 is known, all the coefficients can be determined.

Fig. 3 .

 3 Fig. 3. Schematic representation of assembled augmented matrix of the system (* :

Fig. 4 .

 4 Fig. 4. The variation of resonance frequency with the number of terms K1. The electrical impedances for three different cases -K1 = K2 = 1 (Fig. 5(a)), K1 = 5, K2 = 2 (Fig. 5(b)), and K1 = 21, K2 = 2 (Fig. -are shown and compared with those of the FE method in Fig.5. To measure the discrepancy, a mean absolute percentage deviation (MAPD) metric[START_REF] Wang | Monitoring of the strength gain of concrete using embedded PZT impedance transducer[END_REF] is introduced and defined as

Fig. 5 .

 5 The modulus of electrical impedance of the piezoelectric disk obtained by using different number of terms. (a) K1 = 1, K2 = 1; (b) K1 = 5, K2 = 2; (c) K1 = 21, K2 = 2. Solid lines refer to the results of the FE model and chain dotted lines refer to the results of the analytical model.The mode shapes of total displacement U,

Fig. 6 .Fig. 7 .

 67 The total displacement and electric potential of the piezoelectric ceramic disk at the antiresonance frequency of the first radial (R1) mode. (a)-(b) indicate results of the FE model; (c)-(d) indicate results of analytical model with K1 = 21 and K2 = 2. The total displacement and electric potential of the piezoelectric ceramic disk at the antiresonance frequency of the first thickness (T1) mode. (a)-(b) indicate results of the FE model; (c)-(d) indicate results of analytical model with K1 = 21 and K2 = 2.

Fig. 8 .

 8 Fig. 8. Resonance frequencies vs diameter-to-thickness ratios D/T of PZ27 disks with constant thickness 1.13 mm. Gray points, scaled by normalized amplitude, are results of the analytical model and black circle are results of the FE model.The two horizontal dotted lines, labeled as k33 and kt in Fig.8, are results from the one-dimensional (1D) KLM model[START_REF] Krimholtz | New equivalent circuits for elementary piezoelectric transducers[END_REF]. As D/T ratio approaches 0, the piezoelectric cylinder turns into a long rod with small section poled along its largest dimension, which can be modelled by KLM model with electromechanical coupling coefficient

Fig. 9 .

 9 The modulus of electrical impedance of a long piezoelectric rod (a) and a thin piezoelectric disk (b). Solid lines refer to the results of the present analytical model and chain dotted lines refer to the results of KLM model.

Table 1 ,

 1 where ρ is the density, E c αβ are the elastic

	stiffness constants under constant electric field, i e α are the piezoelectric constants, S ij ε
	are the dielectric constants under constant strain, m δ and e δ are respectively the
	mechanical loss factor and dielectric loss factor.
	Table 1
	Soft PZ27 piezoelectric material properties [36]
	ρ (kg•m -3 )	11 E c (GPa)	12 E c (GPa)	13 E c (GPa)	33 E c (GPa)	44 E c (GPa)	66 E c (GPa)
	7800	147.4	104.9	93.7	112.6	23.0	21.25
	31 e (C•m -2 )	33 e (C•m -2 ) 15 e (C•m -2 )	ε	r	11 S	ε	r	33 S	m δ (%)	e δ (%)
	-3.1	16.0	11.6	913.7	1129.7	1.35	1.7

Table 2

 2 Resonance and antiresonance frequencies of the piezoelectric disk from 0 to 2.4 MHz

	Model						Resonance frequency (MHz)
	FE		0.124 0.315 0.485 0.636 0.762 0.864 0.949 1.018 1.062 1.300
	Analytical		0.126 0.315 0.486 0.634 0.762 0.864 0.949 1.020 1.065 1.300
	(deviation (%))	(1.61) (0.00) (0.21) (0.00) (0.00) (0.00) (0.00) (0.20) (0.28) (0.00)
	FE		1.392 1.481 1.560 1.649 1.736 1.825 1.930 2.051 2.187	̸
	Analytical		1.392 1.482 1.561 1.650 1.736 1.824 1.929 2.050 2.187	̸
	(deviation (%))	(0.00) (0.07) (0.06) (0.06) (0.00) (0.05) (0.05) (0.05) (0.00)	̸
	Model						Antiresonance frequency (MHz)
	FE		0.146 0.324 0.494 0.645 0.772 0.875 0.961 1.032 1.075 1.308
	Analytical		0.147 0.325 0.494 0.645 0.772 0.875 0.961 1.033 1.078 1.308
	(deviation (%))	(0.68) (0.31) (0.00) (0.00) (0.00) (0.00) (0.00) (0.10) (0.28) (0.00)
	FE		1.404 1.512 1.578 1.667 1.804 1.891 1.953 2.070 2.201	̸
	Analytical		1.404 1.512 1.578 1.668 1.804 1.891 1.953 2.068 2.199	̸
	(deviation (%))	(0.00) (0.00) (0.00) (0.06) (0.00) (0.00) (0.00) (0.10) (0.10)	̸
	10	4	R1
	10	3									T1
	10	2
	10	1
			0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2	2.2
								(a)
			0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2	2.2
								(b)

  )

		13	55	33	31	15	33
	1
			15 33	33 11
					2
			33 33	33
	2	2
			15 33	33 11
		55 33	+	15 33		2 ω ρ	33
	3						2
		15 33		33 11	15 33	33 11

  )

	2 3 = n n b b	( 55 11 15 33 ( 15 55 31 15 33 + E E S S e c e e 13 33 11 ) ( ) 55 15 33 15 33 2 15 15 31 13 + + + ( 31 33 11 55 15 + E E S S S E S e e c c e e ε ε ε -+ -+ E c c e e e e e ) ) 15 15 33 e c e e 55 E c e e ε ε ε   = - + -	   	k	n	-	(	) ω ρ 2 15 e e k 15 31 e +	n	(A. 3)
			n	4			(	2 33 33 33 E S -S	)	15 2 n S	15 33 S S	S	2 ω ρ
								15 33	33 11	15 33	33 11
						(	13 E			55 E	)	11 S	15	(	31	+	15	)
			1 n							(	S 15 33	S 33 11	)	n
												55 11 E S	+	2 15
			n	2					2 n	(	S 15 33	S 33 11	)
								33 11	15 33	2 ω ρ	11
				3
								15 33	33 11	15 33	33 1
	13 E c c 55 15 33 E S e ε ε + = 1 n a	11 S -	31 33 11 S e e e 15 15 + n e k ε +	+	31 e	11 15 E c e k +	n
	n a	2	55 31 E c k e e 15 n = -+	+	31 e	2 15 ρω e k +	n
	n a	3	= -	55 11 15 33 E S S c k e 2 n ε ε	2 15 33 11 S e + e ε -
	4 = n a		33 11 15 33 E S S c e ε ε	+ -	15 33 33 11 S e e e ε	+	13 31 E c c 55 15 E e e + +	-	2 ω ρ 15 33 S e e 33 11 11 S S ε ε ε -	k	2 n
	1 n b	= -	(	31 e	) e k 15 11 15 E e c +	n

Acknowledgements

The work was financially supported by the China Scholarship Council (CSC) through the cooperation program UT-INSA (France).

The authors thank Meggitt A/S Denmark for providing the Ferroperm Piezoceramic samples.

Appendix

The constants ami and bmi (i = 1,2,3,4) appearing in Eq. ( 12) are defined by (

) ( ) e e e

and the other coefficients bi, ci and di (i = 1,2,3) are defined as (

)


Florian Lardeux
email: florian.lardeux1@univ-lr.fr

Sylvain Marchand

Petra Gomez-Krämer

Low-complexity arrays of contour signatures for exact shape retrieval

Keywords:

a b s t r a c t

We propose a framework for a fast exact shape retrieval called Low-complexity Arrays of Contour Signatures. The purposes are to match a shape against a database in constant time and to retrieve correct shapes very close to the query, while the latter may have undergone rigid transformations and noise. We present a shape signature based on prior works as well as a compact characterization of such signatures, a system of associative arrays allowing a short search time for retrieval and a technique of pairwise alignment. This method shows a good resilience to perturbations and is performed in a constant computational time.

Introduction

Shape recognition, classification, retrieval or even just building tools for shape analysis can still be a challenge nowadays, as the exigencies of technological effort s continuously increase. There are many ways to formulate goals within this field: classifying shapes visually resembling each other, finding matching parts of shapes, perform a matching that is invariant or robust to transformations, etc. A significant part of the literature focuses on shape classification considering occlusions, perturbations due to noise, scaling, rotation and even sometimes large deformations. Depending on the application, the intra-class variance in a dataset is sometimes such that it can be considered relevant to admit that some shapes in the same class are not the same. Therefore a large part of the existing state of the art cannot discriminate visually similar but differently labeled shapes, as they admit large amounts of deformations. Also, most of the current literature revolves around the use of direct descriptor comparisons, which can lead to high computational times for retrieval.

In this paper, we decide to tackle the problem of exact shape retrieval in a large database, in constant computational time.

For instance, these shapes can be contours extracted from quasi-flat objects such as ancient coins [START_REF] Lardeux | Multi-light energy map[END_REF] . Ancient coin databases can be very large, i.e. the Online Coins of the Roman Empire (OCRE) database 1 totals over 10 0,0 0 0 specimens. The objective is to retrieve only coins that are produced by the same matrix. This is a problem of exact shape retrieval different from typical shape retrieval approaches, which would retrieve all coins of the same type.

Hence, two shapes are considered here to be the same if one can be mapped onto the other by applying a rigid transformation or noise, as presented in Fig. 1 .

The intuition for this work was given by a conjunction of ideas from both the audio field and the shape recognition community. The work of [START_REF] Cui | Curve matching for open 2D curves[END_REF] is a very simple and compact way of describing the shape by its curvature. However, the cross-correlation is not robust and quite slow. Wang et al. introduced in [START_REF] Wang | An industrial strength audio search algorithm[END_REF] the Shazam system for musical audio recognition, where retrieval is performed using a hash map to efficiently store feature pairs of the signal as well as a binning technique to ensure the alignment between similar parts. We therefore investigate this idea in the context of shape recognition, without the need for a hash key for the creation of the feature pairs. Dealing with small deformations involves nonetheless the use of a robust signature. Since the cumulative sum [START_REF] Cui | Curve matching for open 2D curves[END_REF] accumulates noise along the signature, this representation is not at all robust to perturbations. Our signature is a modified version of the Triangle Centroid Distances (TCD) [START_REF] Yang | A novel method for 2D nonrigid partial shape matching[END_REF] .

This article involves the following contributions: 1) A signature based on a modified version of the TCD [START_REF] Yang | A novel method for 2D nonrigid partial shape matching[END_REF] , which makes a simple yet very robust 1D representation of the shape. 2) A selection of feature points and feature pairs extracted from the signature. 3) An efficient storage system via associative arrays. 4) A general algorithm of constant time retrieval using pairwise geometric alignment. We call this whole method Low-complexity Arrays of Contour Signatures, shortened in LACS.

The use of associative arrays for exact shape matching resembles geometric hashing. However, geometric hashing was originally Fig. 1. Exact shape retrieval. A shape is retrieved if and only if there is a rigid transformation or noise between it and the query. If the shape seems visually to belong to the same class, but is in fact geometrically different, it will not be retrieved. developed in computer vision for matching point features [START_REF] Wolfson | Geometric hashing: an overview[END_REF] , such as vertices of shapes, and not shape signatures.

Related works on shape recognition are listed in Section 2 . The definition of the signature is presented in Section 3 . In Section 4 , the feature points and the feature pairs are described and the process of storage to form a database is described in Section 5 . Section 6 explains the mechanism of pairwise geometric alignment to ensure geometric consistency in the retrieval. Finally, experiments on robustness and runtime are presented in Section 7 . A conclusion and discussion of our work are outlined in Section 8 .

Related work

There are two main perspectives to extract information: region based and contour based. Region-based methods, such as Fourier descriptors [START_REF] Zhang | Shape-based image retrieval using generic Fourier descriptor[END_REF] or image moments [START_REF] Liao | On image analysis by moments[END_REF] , provide a global shape description by computing values related to the shape as a whole. However, these techniques suffer from the global descriptive data since they can analyze only parts of the shape.

Contour-based methods are widely considered in the literature, because they are invariant to lighting conditions and do not take into account colors or textures [START_REF] Riemenschneider | Using partial edge contour matches for efficient object category localization[END_REF] . The main strategy to characterize the contour is the use of a shape descriptor, which can characterize globally or locally the contour. The corner stone of shape similarity techniques is the shape context descriptor [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF] , which locally features the distribution of contour points by the means of radial histograms. Although it allows to deal with scaling, rotation and deformations, the similarity cost is rather heavy and cannot cope with partial occlusions. Many works are based on this technique and propose the use of Dynamic Programming to match points [START_REF] Daliri | Robust symbolic representation for shape recognition and retrieval[END_REF] , a hierarchical strategy to deal with affine and projective transformations [START_REF] Jia | Hierarchical projective invariant contexts for shape recognition[END_REF] or the use of the earth mover's distance to perform comparisons [START_REF] Shu | A novel contour descriptor for 2D shape matching and its application to image retrieval[END_REF] . Chang et al. [START_REF] Chang | Partial shape matching and retrieval under occlusion and noise[END_REF] propose a descriptor for open and closed shapes based on the Triangle-Area representation. Their work efficiently deal with scaling by combining several scales and can tackle occlusions. Similar properties can be found in [START_REF] Yang | Invariant multi-scale descriptor for shape representation, matching and retrieval[END_REF][START_REF] Kun | Shape matching based on multi-scale invariant features[END_REF] , which propose multiscale descriptors based on invariant aspects of shapes. They also demonstrate invariance to articulated deformation and intra-class variations. Wang et al. [START_REF] Wang | Shape matching and classification using height functions[END_REF] introduce another type of descriptor relying on height functions. The works of [START_REF] Donoser | Efficient partial shape matching of outer contours[END_REF] and [START_REF] Riemenschneider | Using partial edge contour matches for efficient object category localization[END_REF] use a descriptor which can efficiently match parts of contours. The general method uses a fixed sampling of the contour from which they extract angles between chords within the shape. Yang et al. introduced the Triangle Cendroid Distance (TCD) descriptor [START_REF] Yang | A novel method for 2D nonrigid partial shape matching[END_REF] , based on distances between each point of the shape and Fourier descriptors. They achieve fine results in both part-to-whole and part-to-part matching while dealing with significant non-rigid deformations. However, their matching technique is not invariant to the position of the starting point and therefore needs as many comparisons per shape as there are points in it.

The work of [START_REF] Zheng | Fourier transform to group feature on generated coarser contours for fast 2D shape matching[END_REF] also propose the use of Fourier descriptors with a multiscale technique using group features, which uses combination of spatial features to improve matching accuracy. Finally, to refine shape description, the work of [START_REF] Shen | Shape recognition by bag of skeleton-associated contour parts[END_REF] combines shape contour and skeleton by analyzing the geometry of contour parts in relation with their associated skeleton.

Alternatively, several other works use signature-based representations to describe the shape, which synthesize the contour information into a compact representation. Commonly used signatures are, for instance, centroid distances, the chord length signature, the area function, the cumulative angular signature and the turning angle signature. Curvature is a paramount element in human shape recognition [START_REF] Attneave | Some informational aspects of visual perception[END_REF] . Cui et al. [START_REF] Cui | Curve matching for open 2D curves[END_REF] formulate a scale-invariant signature based on this principle. The abscissa is calculated as the integral of unsigned curvatures. This mathematically shows to be indeed scale-invariant in theory but in practice fails to take into account the changes in sampling of the contour when perturbations are applied to the image itself. To compare two shapes, they use normalized cross-correlation. In [START_REF] Kaothanthong | Distance interior ratio: a new shape signature for 2D shape retrieval[END_REF] , a Distance Interior Ratio signature is introduced. It is invariant to rotation and scaling and yield promising retrieval results on popular datasets.

Various frameworks can be formulated for the matching step, depending on what the desired goal is. Research works describing novel methods for shape representation generally adopt distances or similarity measures such as the edit distance [START_REF] Daliri | Robust symbolic representation for shape recognition and retrieval[END_REF] , the earth mover's distance [START_REF] Grauman | Fast contour matching using approximate earth movers distance[END_REF] , the χ 2 distance [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF] , or dissimilarity costs [START_REF] Yang | A novel method for 2D nonrigid partial shape matching[END_REF][START_REF] Donoser | Efficient partial shape matching of outer contours[END_REF] . In any of the aforementioned cases, pairwise comparisons must be performed if one desires to find the best match among a database of shapes which leads to a linear complexity when performing retrieval. Several authors propose alternatives to standard comparisons that can involve space search reduction via statistical tools [START_REF] Lim | Fast shape matching using statistical features of shape contexts[END_REF] , or the use of Bayesian inference [START_REF] Zhang | Bayesian clustering of shapes of curves[END_REF] . However, the speed gain is achieved for each comparison, and not for the database as a whole, which is still a limiting factor for a fast retrieval. Other works rely on learning to perform classification. The classic method of Bag-of-Words model is used by [START_REF] Wang | Bag of contour fragments for robust shape classification[END_REF][START_REF] Zeng | Curvature bag of words model for shape recognition[END_REF][START_REF] Shen | Bag of shape features with a learned pooling function for shape recognition[END_REF] . Bai et al. [START_REF] Bai | Shape vocabulary: a robust and efficient shape representation for shape matching[END_REF] extended this principle by using both local information and global properties of the shape to gain discriminative power for the recognition. The same idea of contour segments is investigated in other works such as [START_REF] Sun | Classification of contour shapes using class segment sets[END_REF] which uses Bayesian inference for learning. Learning methodologies can be very effective for shape classification, but require a specific training set, which cannot account for large databases of very different shapes.

The shape as a signature

Our choice is to convert the shape into a 1D signature so as to get a compact representation which contains relevant information of the shape. Attneave [START_REF] Attneave | Some informational aspects of visual perception[END_REF] showed in 1954 that the human mind seems to privilege a very small, but really precise quantity of data, gathered in its boundary curvature.

In order to comply with this psychological background, a natural choice is to study signatures which feature this kind of attributes. Most naive signatures [START_REF] Zhang | A comparative study of Fourier descriptors for shape representation and retrieval[END_REF] are not invariant to scale and features such as curvature also suffer heavily from noise [START_REF] Cui | Curve matching for open 2D curves[END_REF] . Namely, the main issue with the computation of a good signature is its stability under perturbations.

Description of the signature

We modify the TCD descriptor designed by [START_REF] Yang | A novel method for 2D nonrigid partial shape matching[END_REF] and extract a signature from it.

The TCD is based on the Fourier transform of distances between points of the shape. Namely, if the shape has n points, for each pair of points P i and P j , j = i , let g i j be the barycenter between P i , P j and the total center of gravity of the shape G . Then, the distance Fig. 2. Graphical explanation of the TCD. d i j is the distance between P i (red) and g i j (green). G is the center of gravity of the whole shape. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) d i j is defined as the distance between P i and g i j . Fig. 2 provides a illustration of this calculation.

Each distance is normalized by the maximum of each row, giving Eq. ( 1) .

∀ (i, j) ∈ [0 , n -1] × [0 , n -2] , ˆ d i j = d i j max i d i j (1)
Finally, a Fourier transform is applied for each column of (d i j )

for 16 frequencies, from which the absolute value is taken, which leads to the final descriptor

( f w j ) ∈ M 16 ,n -1 ( R + ) in which each el-
ement is given by Eq. ( 2) .

∀ j = 0 , . . . , n -2 , ∀ ω = 0 , . . . , 15
f ω j = 1 (N -1) • max n d n j • n = j d n j • exp - 2 iπ nω N -1 (2)
where i 2 = -1 and | . | is the absolute value.

The difference in description compared to the original TCD is that we consider the d n j to be the direct norm of the vector --→ P n P j . This formulation is much simpler and yields better results because it deals with exact positioning. Our principal signature, which is sufficient for a large number of cases, is the case of zero frequency in Eq. ( 2) , and is simply the ratio of the average and the maximum of the d n j .

f j = 1 (N -1) • max n d n j • n = j d n j = d n j max n d n j ( 3
)
This signature is robust to rigid transformations as well as noise and characterizes the shape. As explained below in Section 5.3 , one can nonetheless use at least the first few rows of ( f w j ) to improve the robustness. It should be noted that each shape is read forward and backward, since there is a priori no way to select a stable direction. The query shape will thus have two signatures to compare.

Behavior of the signature

This simple formulation is very efficient as a signature since it is, in theory, invariant to shift, scale, rotation, starting point (resulting in a shift in abscissa) and very robust to noise, being the low frequency of the distance signal. Here we focus on the formulation presented in Eq. ( 3) . Three types of behavior regarding transformations are cited below: adding noise the to shape, scaling the source and rotating the source (changing the starting point). For illustration, the first shape of the Kimia99 dataset [START_REF] Sebastian | Recognition of shapes by editing shock graphs[END_REF] , displayed in Fig. 4 , is used. All original images have an average resolution of around 120 × 120 . A binary blob is extracted from the shape and the contour is extracted as a contiguous sequence of pixels.

Scaling the image

Scaling is applied to the image prior to the shape extraction and the calculation of its signature. The larger the image, the more precise the shape details. Fig. 3 presents the two versions of the signature of the same shape having undergone a scaling by 2. It is paramount that the signature is as precise as possible so that its extrema stay very stable (see Section 4 ). The work in [START_REF] Yang | A novel method for 2D nonrigid partial shape matching[END_REF] uses a uniform interpolation of 100 points whereas we choose 10 0 0 to gain precision.

Rotating the image

Rotation changes the location of the starting point within a closed shape. Given a binary blob image from which the contour is extracted, there is no way to determine robustly which point to start the contour with, and therefore it changes with rotation. This is particularly visible in Fig. 4 , where the red dot denotes the starting point of the shape. Fig. 5 reveals that the signature shifts accordingly. The only major issue with this phenomenon is that it can hide a potentially relevant feature point being located at the break point that would appear as a peak or a trough in the signature.

Adding noise to the shape

In order to simulate noise, we simply add a Gaussian noise to both coordinates of the shape pixels independently. Fig. 6 shows the impact of adding noise on the signature. A Gaussian noise of standard deviation 1.0 is rather strong for the size of this shape, therefore it can be considered as a limit case. The signature can be very corrupted by noise, the low frequency signal is used. This is considered to be a limit case, since we do not intend to deal with large deformations.

Feature points and feature pairs

To reduce both spatial and time complexities, feature points are extracted from the signature values and organized in pairs. That way, we ensure to keep the model very compact while maintaining a geometrical consistency.

Feature point selection

Feature points of our signature are defined as the points which exhibit a distinct variation in signature.

The feature points are chosen to be the local extrema of the signature of coordinates (x j , f j ) such that:

P feature = (x j , f j ) | f j > f j-1 and f j > f j+1 or f j < f j-1 and f j < f j+1 (4)
Fig. 7 shows these feature points in our example in both the signature and the shape. A simple analysis of Eq. ( 3) reveals that extrema in the signature appear where points are the closest and the furthest from the rest of the shape. This may occur at various positions in the shape, but this set of feature points noticeably include high curvature points. However, they carry more information than the sole extrema of the curvature. This is at the heart of a relevant and compact representation of the shape for an efficient recognition; shapes that carry very little curvature extrema cannot be analyzed properly. Here, on the contrary, there is a nice distribution of peaks and troughs that enable a rich feature representation of the shape. Finally, by keeping only the extrema, a large part of the information carried by the other parts of the signature are eliminated. This pruning is relevant because we postulate that the informative content is gathered within the relative positions of these extrema.

Generating pairs of feature points

Out of the selected feature points, our idea is to generate feature pairs. The main argument behind such a choice is the following: points by themselves do not carry sufficient information; there is no way to compare two points if there is a shift between them. Pairs of points add structural information about the signature and the shape. Because the signature f is very robust to rigid deformations, the pairs will remain very stable.

The generation of pairs is done by reading the feature points in the signature from left to right. For each point, the algorithm looks ahead within a window of size m inside which the second point is selected. An example of such pairs can be found in Fig. 8 . In practice, we choose m = 20 .

To handle the case of the rotation of a closed shape and thus a shift of the signature, one needs to take into account the points in a circular manner. The list of n feature points is appended with the first m elements, making it a list of n + m feature points. This way, when processing the last points of the list, we create pairs with the first feature points of the list, thus creating a continuity in the generation of pairs. In the same fashion, we calculate another list of feature pairs for the backward signature.

Storing data in associative arrays

We choose to use the pairs to store information in a system of associative arrays that serves as a tool for an efficient retrieval method (see Section 6 ). The main idea behind this form of storage is that one can store a whole database of feature pairs in it, each belonging to a certain shape. In essence, this is very close to the concept of hash tables. Some works have used hashing for shape recognition and clustering [START_REF] Lamdan | Geometric hashing: A general and efficient model-based recognition scheme[END_REF][START_REF] Kaplan | Locality sensitive hashing for efficient similar polygon retrieval[END_REF] , or for finding unusual shapes [START_REF] Wei | Efficiently finding unusual shapes in large image databases[END_REF] . However, our technique allows to retrieve shapes among a large database and does not require learning any clustering pattern. Although the hashing technique has also been used in the Shazam algorithm [START_REF] Wang | An industrial strength audio search algorithm[END_REF] by which our method was inspired, the difference is that we do not explicitly create hash functions per se , because the associative functions described below are not injective.

Associative array functions

An associative array is a set of buckets indexed by a key function and containing a certain value. Let h = { h 0 , h 1 , . . . , h p-1 } (p ∈ N * ) be the set of associative arrays. Each array does not have a fixed number of buckets. In the following, each array is assimilated to its key function so that

∀ i = 0 , . . . , p -2 , h i ∈ Z R .
Every time a pair is stored in h , the bucket indexed by the key might not exist and is therefore created in the process. In prac-Fig. 7. Feature points in our example shape. Note that these extrema are closely related to the high curvature points, but give however far more information than a simple curvature signature.

d i = h i (d i ) 4 if i = p -1 then 5 if L [ ˆ d i ] does not exist then 6 L [ ˆ d i ] ← empty list 7 L [ ˆ d i ] ← ---- append (x i , idx ) 8 else 9 if L [ ˆ d i ] does not exist then 10 L [ ˆ d i ] ← empty array 11 L ← L [ ˆ d i ]
tice, we assign an index to the shape ( idx in Algorithm 1 ), which can be of any form as long as it allows to identify the shape.

For each pair, there is a descriptor

D = { d 0 , d 1 , • • • , d p-
h i (d i ) = a i • d i , ∀ i = 0 , . . . , p -1 .
Because the signature is deemed sufficiently stable, one can discriminate the ordinates easily, hence the choice of the first two components. Moreover, the third component of D allows to neglect the impact of shift along the signature abscissa, thereby making the descriptor invariant to the starting point. The descriptors are thus very reliable, provided the signature f is sufficiently robust.

Structure

As evoked, h is a system of associative arrays. Each bucket (not empty) of the first array h 0 contains one independent array, which is indexed by the second key function h 1 , and so on. This system of nested arrays is akin to a dense tree structure, in which there the retrieval is progressively refined. Fig. 9 explains graphically the structure of our system and presents a query example. In our case, we use three different keys and therefore the system has a depth of 3.

The last array buckets should contain the identification of the shape. We store a list of tuples consisting of the abscissa of the first point x 1 and the index of the shape idx . Therefore, each time a new entry is added to the whole array ( Algorithm 1 ), it is possible to append its tuple (x 1 , idx ) to a list already containing previously added tuples. The associative array h i does not handle collisions and the whole storage system allows to retrieve several potential shape indices for each query pair ( Algorithm 3 ). Also,

( f 1 , f 2 , x 1 -x 2 ) ,
the result is found by successively finding the corresponding bin of the arrays, as long as it exists. The last array contains the result of the retrieval, which is the abscissa of the first point of the retrieved feature pair x 1 initially stored ( Algorithm 1 ) and the index of the corresponding shape idx .

Fig. 10. Example of H

= { h u v | u ∈ [0 , 2] v ∈ [0 , 2]
} . Each h u is a system of associative arrays composed of v nested arrays. The retrieve values (x 1 , idx ) may not concur, as they are pulled from a different signature (even though they characterize the same shape, they are independently processed).

if no corresponding pair was found (which should often occur if the shape is unknown in the database), then the value None is returned.

Parallel systems of arrays

It is possible to use several systems of associative arrays at once. Namely, we can use several signatures, all extracted from the descriptor of [START_REF] Yang | A novel method for 2D nonrigid partial shape matching[END_REF] and processed separately, all in the same fashion, both for storage and retrieval. Let H be a system of associative arrays. The process of retrieval is to consider a set of signatures characterizing the shape at several levels of description f j (ω = 0) , f j (ω = 1) , . . . and perform retrieval on each pair of each signature. The complexity of retrieval increases linearly with the size of H, but adding more systems of arrays leads to very robust results. An example of such set of arrays is presented in Fig. 10 . In the following, we will continue to use only one array to explain the technique to avoid heavy notations.

Retrieval

This section explains the complete mechanism of retrieval. It consists of two general steps: finding the matches in the database and aligning them with a binning technique.

Pulling correspondences from h

The simple scenario to consider here is having a query Q , which is a collection of pair descriptors for a query shape, and compar- ing it with the shapes already processed and stored in h . Let S be the database of shapes stored in h . The shape k in the database is

S(k ) , k ∈ { 0 , 1 , . . . , n -1 }
, where n is the number of shapes in the database.

The algorithm goes through all pair descriptors in Q , and for each of them, we retrieve a list of tuples of abscissa and index that we write { ( x, k ) } for clarity. If the retrieved value is None, the next descriptor of Q is processed, ending up retrieving all the potential shape candidates for each investigated pair. The retrieved values lack information (only one abscissa value x along with its corresponding shape index k ) and may be numerous for each query pair. Therefore an additional strategy is employed to select the relevant ones and prune the acquired set.

Pairwise geometric alignment

If x and x Q are respectively the abscissa in the current retrieved tuple and the first point of the current pair in the query, we define:

δx x -x Q (5)
Eq. ( 5) refers to the shift needed to map x Q onto x . An example of such calculation is provided in Fig. 11 . For any part of the signature which matches the query, the resulting values of δx should be close to another one. The values of δx are quantized ( Eq. ( 6) ) and we calculate histograms of δx per shape of index k in the database. The bins of the histograms range from δ min = -10 0 0 to δ max = +10 0 0 and have N bins bins. The resulting binning is expressed by Eq. ( 6) . The number of bins is again a parameter for the finesse of the matching; values ranging from 100 to 500 give fairly good results and here we choose 200 bins. We determine the quantized values of δx on the fly for each retrieved tuple, and we successively identify the score of the current shape score [ k ] as the value of the largest bin. The scores are stored in a dictionary indexed by the shape indices.

x = N bins • δxδ min δ max -δ min [START_REF] Zhang | Shape-based image retrieval using generic Fourier descriptor[END_REF] Fig. 12 illustrates the behavior of such histograms. If the original shape is matched against itself, the correspondences are perfect, and therefore the retrieval yields a maximum number of pairs. This number is far larger than the number of pairs in the query because it is likely to find several results corresponding to one input descriptor. As explained in Section 3.2 , noisy shapes disturb the signature signal and can thus corrupt the retrieval and rotation is the reason for an important shift in the signatures which can be seen in the histogram. There is a significant drop in the number of pairs retrieved, but nonetheless a very coherent shift is found.

Selection of the best shape

The aforementioned pipeline of techniques is compiled into a single Algorithm 3 . There are essentially two passes to compute the scores: one for the forward query and one for the backward query. At each iteration, we retrieve the correspondences from the database and update the score and the best shape associated using the steps described in Section 6.2 . We only keep the largest score for each bin of each shape. The best score, corresponding to the largest number of pairs retrieved over all shapes in S is

n best = max k (score [ k ]) .
Finally, the index of the best shape is given by k best = arg max k (score [ k ]) .

Computational complexity

The descriptors and the quantization produce collisions during storage. This allows to get more than one possible shape corresponding to the query and also several pairs belonging to the same shape. Of course, the main issue is that it causes a linear computational complexity instead of a purely constant one. In theory, given a set of associative arrays, each having b i buckets, the mean number of collisions can be assumed to be uniformly distributed and equals to:

γ = N/ p-1 i =0 b i ( 7
)
where N is the number of pairs stored. The number of pairs N is substituted by the total number of shapes in the database n since the former increases linearly with n . Using Eq. ( 7) , the complexity becomes:

C(n ) = O(γ • N Q ) = O n b p • N Q (8)
where N Q is the number of pairs in the query and if we suppose i b i ≈ b p . Eq. ( 8) shows that there is a linear relation between the complexity and the number of pairs N, but there is also a factor which, if well-chosen, can severely diminish the overall computational time. The value of b depends on the descriptor definition and also on the parameters A ; a carefully chosen set of parameters thus allows to limit collisions, to ensure that γ is relatively constant when N (and the number of shapes) varies and to only pick pairs closely related to the ones investigated as queries. We will demonstrate experimentally in Section 7.2 that we manage to keep the runtime low and γ quasi-constant.

Experiments and results

To assess the performance of exact shape retrieval, studies on both robustness to perturbations and retrieval time are provided. Our technique is compared to several methods from the state of the art, including [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF] as a standard baseline strategy with descriptors, [START_REF] Cui | Curve matching for open 2D curves[END_REF] to evaluate the use of cross-correlation and [START_REF] Yang | A novel method for 2D nonrigid partial shape matching[END_REF] because their method inspired our signature. In the graphs described below, our method is labeled "LACS u ", where u corresponds to the number of parallel systems of arrays used (see Section 5.3 ).

Robustness

This subsection presents an analysis of the resilience with respect to scaling, rotation and noise. These perturbations were generated by a process explained in Section 3.2 . As a general protocol for this study, each transformed shape query is matched against the original shapes. Parameters for the perturbations are presented in Table 1 . We choose to use the Kimia99 dataset [START_REF] Sebastian | Recognition of shapes by editing shock graphs[END_REF] as a database for two reasons: it features a reasonable variety of inter-class shape complexities, each class sometimes containing very similar shapes, albeit without being a too large dataset, which leads to prohibitive computational times, especially for [START_REF] Cui | Curve matching for open 2D curves[END_REF][START_REF] Yang | A novel method for 2D nonrigid partial shape matching[END_REF][START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF] . Fig. 13 shows some original and perturbed shapes of the Kimia99 dataset. The final values are accuracy scores corresponding to the bull's eye score for the retrieval of each independent shape, the retrieval per class and per parameter of perturbation. Results in Table 2 show that our method outperforms both [START_REF] Cui | Curve matching for open 2D curves[END_REF] and [START_REF] Yang | A novel method for 2D nonrigid partial shape matching[END_REF] , even with one associative array. Considerable improvements in robustness can be made by adding more associative arrays, but at the cost of more computational time. Shape context [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF] achieves impressive results of robustness to noise, but is much less resilient for scaling and rotation. The TCD is designed to deal with non-rigid deformations, which obviously hinder its discrimination for very close shapes. The accuracy per class is however still very high in general, meaning that mistakes were mostly made between very close shapes ( Fig. 14 ). Cross-correlation of [START_REF] Cui | Curve matching for open 2D curves[END_REF] suffers from a very poor stability, as predicted. Our method achieves very good outcomes in general for the accuracies per transformation parameter. Scaling mistakes correlate with factors higher than 1.0. Image rotation involves drops in robustness that are more likely to appear for angles θ ≡ 0 ( mod π / 2) . Finally, there is a rapid decline in robustness for higher noise standard deviations for all methods, except for the shape context, but we still manage to reach 68% for σ = 1 . 0 .

We show some further results in Table 2 for the MPEG-7 dataset [START_REF]The Moving Picture Experts Group[END_REF] and the Tari-10 0 0 dataset [START_REF] Asian | An axis-based representation for recognition[END_REF] . They are usual datasets for the evaluation of shape retrieval methods and contain respectively 1400 and 1000 shapes. The results for LACS and TCD are comparable to those obtained on the Kimia99 dataset. As shown below, shape context and the method of Cui et al. have a quadratic time complexity. Hence, it is not possible to execute them on such a large dataset in a reasonable time interval.

Runtime analysis

A runtime analysis is performed to evaluate the behavior of the associative arrays when changing the size of the database and the size of the query. We use the MPEG-7 dataset [START_REF]The Moving Picture Experts Group[END_REF] , which contains 1400 shapes, because we want to evaluate the runtime for large databases. For the first experiment, we select a specific query containing N Q = 380 pairs, a value close to the mean number of pairs N Q = 374 . We perform successive retrievals against a database in which we iteratively increase the number of shapes n ; the same process is followed for every method. The second experiment is done by using the largest database (1400 shapes are stored) and by computing the retrieval time for an increasing number of pairs in the shape. To do so, we take the shape with the largest number of pairs and perform a process similar to the first experiment by adding pairs. In each experiment, we do not take into account the creation of the descriptors, the signatures or the associative arrays; the database is already already built and we examine only the retrieval time. Computations were run with an Intel® Core TM i7-7820HQ CPU @ 2.90GHz with Python 3.6.5 and the results are presented in Fig. 15 .

Although we indeed observe a linear tendency, the curve is very flat, making the computational time of technique akin to a constant complexity. When the database is very small, there is a high probability of retrieving nothing (None in Algorithm 2 ) and so the retrieval time is negligible. This explains the soaring at first, which evens out at around 100 shapes. The average number of collisions (calculated with Eq. ( 7) ) remains between γ = 1 . 31 and γ = 2 . 92 , with a general mean of γ = 2 . 27 . This implies the necessity to check, in average, 2.27 values r in Algorithm 3 . The result is a very low retrieval time which spans between 3 ms and 13 ms, with an average of 10 ms. One can extrapolate that only 1 s is required to perform a retrieval with a query of 380 pairs among a database of around n (t = 1 s ) = 350 0 0 0 shapes. Table 3 demonstrates the speed of our method with respect to more traditional approaches. In av-Algorithm 2: Request on the associative array.

Input

: h : associative array erage, it takes less than 1 s to perform retrieval for any number of parallel associative arrays. By comparison with direct matching, our method is 10 to 10 0 0 times faster for the range of n analyzed. We also present the number of collisions with respect to the n , which exhibit a linear correspondence between γ and the computational time. Fig. 16 features the evolution of runtime against N Q as well, and explicitly shows a linear tendency. The number of pairs N Q is positively correlated to the complexity of the shape, as defined by the number of high curvature points and its length.

Since we expect an average of N Q = 374 pairs, the average computational time is of around 100 ms for n = 1400 .

Conclusion and discussion

We introduced the Low-complexity Arrays of Contour Signatures (LACS) system for shape retrieval. The scope was to deal with exact shape recognition in a constant time, regardless of the size of the database. Our method formulates a very robust representation even in the presence of rigid perturbations and noise. Retrieval is performed with the use of a system of associative arrays and simultaneous pairwise geometric alignment, which yields extremely short computational times on large databases compared to usual works on the subject. Our technique only tackles exact matching and so far cannot deal with occlusions, which comes from the design of the signature. The pipeline for retrieval is largely interchangeable, as one may use other shape representations in our system. Further improvements could lead to very satisfactory results in terms of recognition, classification as well as part-to-part matching in a versatile and very fast way.

Fig. 3 .

 3 Fig. 3. Signatures obtained with two different scalings of the shape. There is a slight difference, but overall the signatures are almost identical.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Left and middle: Impact of rotation on the starting point (red dot). The shape here is rotated by π 4 . Right: Shape with Gaussian noise of standard deviation 0.5. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .

 6 Fig. 6. Impact of noise on the signature. Noise in the shape translates in noise in the signature. Strong noise can significantly corrupt the signal and disturb the local extrema locations.

Fig. 8 .

 8 Fig. 8. Examples of pairs of feature points (10 are presented). Pairs add structural information to the positional information of the feature points.

Algorithm 1 :

 1 Adding a pair to associative arrays. Input : h : associative arrays { a 0 , a 1 , . . . , a p-1 } quantization parameters { d 0 , d 1 , . . . , d p-1 } : input descriptor { x 0 , x 1 , . . . , x p-1 } : list of abscissas of the first point of the query pairs idx : index of the shape 1 L ← h 0 2 for i from 0 to p -1 do 3 ˆ

  1 } with p components. Each key function takes as an input respectively each component of the descriptor D. Let A = { a 0 , a 1 , . . . , a p-1 } be a set of p positive values. These values are quantization parameters of the key functions defined below. They virtually allow to increase or decrease the "size" of each bucket in the arrays. This way, one can manipulate those parameters in order to get a finer or coarser storage. Small perturbations can be taken into account by setting a sufficient coarseness. For our application and for any pair of extrema { (x 1 , f 1 ) , (x 2 , f 2 ) } defined in Eq. (4) , we use a descriptor of three components defined by D = { f 1 , f 2 , x 2x 1 } and a set of parameters { a 0 , a 1 , a 2 } such that a 1 = 0 . 1 , a 2 = 0 . 1 and a 3 = 0 . 5 . The values of the arrays are fetched by finding bins corresponding to

Fig. 9 .

 9 Fig. 9. Structure of our system of associative arrays. Each bin of each array, indexed by its key function, contains a sub-array. For a query pair descriptor ( f 1 , f 2 , x 1 -x 2 ) ,

Fig. 11 .

 11 Fig. 11. Example of the calculation of δx for a rotation of π 6 . The given points P 1 , P 2 and P 3 have been shifted due to rotation. Here δx = 166 in each case. Even if there were slight variations (bounded by +2.0 (excluded) with a 3 = 0 . 5 ), those values would still belong to same bin in the histogram.

Fig. 12 .

 12 Fig. 12.Results of the alignment for different perturbations. The first row displays the contour extracted from the transformed shape, the second row shows the corresponding signatures and the third row presents the resulting histograms. It should be mentioned that the noise in this case is quite large and represents a worst-case scenario.

Table 1 Fig. 13 .

 113 Fig. 13. Original and perturbed shapes of the Kimia99 dataset.

{ a 0 Output : list of tuples 1 L 3 ˆ 5 return None 6 L

 01356 , a 1 , . . . , a p-1 } quantization parameters { d 0 , d 1 , . . . , d p-1 } : input descriptor ← h 0 2 for i from 0 to p -1 dod i = h i (d i ) 4 if L [ ˆ d i ] does not exist then ← L [ ˆ d i ]

Fig. 14 .

 14 Fig. 14. Detailed results of the robustness experiments. The left column shows accuracy against the perturbation parameters. The right column shows accuracy per class for each perturbation.

Fig. 15 .

 15 Fig. 15. Left: Retrieval time evaluations with respect to the size of the database for our technique with one array, using same query of N Q = 380 pairs. A linear regression is shown beginning with n = 100 . Right: Retrieval time per collisions γ .

Fig. 16 .Algorithm 3 :do 6 R 7 if R is not None then 8 for all r in R do 9 x = r(0) 10 idx = r( 1 ) 11 δx 12 x 14 bin + = 1 15 if bin ≥ score then 16 k = idx 17 score = bin 18 return k , score

 1636789101111214151617 Fig. 16. Retrieval time evaluations with respect to the size of the query, using the largest database. Runtimes are very low, showing a near constant complexity.

Table 2

 2 Total accuracy of robustness ( % ) by method including various amounts of parallel associative arrays.

		Shape Context [9]	Cui et al. [2]	TCD [4]	LACS (Ours) (number of parallel arrays)
	Kimia99				1	2	3	4	5
	Scaling	66.6	13.7	71.5	80.2	90.3	90.5	92.4	92.8
	Rotation	35.0	13.9	88.9	88.7	94.0	95.0	95.6	96.3
	Noise	95.8	26.8	63.5	66.4	77.3	82.4	85.5	86.9
	MPEG-7				1	2	3	4	5
	Scaling			70.0	86.7	89.9	90.5	91.3	91.1
	Rotation			74.6	83.5	87.4	87.9	88.4	88.6
	Noise			63.7	68.2	76.2	79.1	80.8	82.8
	Tari-1000				1	2	3	4	5
	Scaling			72.6	93.5	95.7	95.8	96.2	96.3
	Rotation			87.4	86.0	89.9	90.4	91.2	91.5
	Noise			75.3	65.6	77.1	80.7	83.0	85.1

Table 3

 3 Runtimes of retrieval (in seconds) for several methods and different amounts n of shapes in the database. Collisions γ for each n is also displayed between parentheses.

	Database size n	Shape Context [9]	Cui et al. [2]	TCD [4]	LACS 1	LACS 2	LACS 3	LACS 4	LACS 5
	1	0.93	0.066	0.010	0.003 (1.31)	0.005 (1.43)	0.011 (1.50)	0.020 (1.53)	0.018 (1.55)
	10	8.95	0.47	0.070	0.002 (1.47)	0.006 (1.50)	0.011 (1.63)	0.022 (1.66)	0.025 (1.68)
	100	86.34	5.19	0.61	0.013 (1.52)	0.009 (1.52)	0.030 (1.61)	0.045 (1.93)	0.121 (2.31)
	1000	823.61	65.82	6.44	0.011 (3.42)	0.011 (3.42)	0.193 (3.91)	0.241 (5.18)	0.901 (6.57)

return L

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.