G Delgado

M Hamdaoui

Topology optimization of frequency dependent viscoelastic structures via a level-set method

Keywords: Topology optimization, Level-set method, Viscoelastic damping, Non-linear eigenvalue problem

Viscoelastic materials follow a liquid-like elastic behavior whose characteristics depend on the excitation frequency. Nowadays, this type of materials represent a high opportunity for vibration damping treatments in the automotive and aeronautic industries, for instance. This work is devoted to the application of the level-set method for topology optimization of viscoelastic structures. We look for the best distribution of viscoelastic material within a reference domain for the design of purely viscoelastic 3D damping structures and 2D viscoelastic damping treatments. In both cases one desires to maximize the structure capacity to dissipate energy measured here by the modal loss factor of the first vibration mode.

Introduction

Viscoelastic damping material behavior occurs in a wide variety of materials and can be characterized by liquid-like elastic behavior. Materials that experience viscoelastic behavior include acrylics, rubber, and adhesives. The characteristics of viscoelastic materials depend on temperature and frequency. Viscoelastic materials have the property of absorbing vibrational energy which makes them very interesting for structural damping applications. The vibrations can be caused by noises that radiate from a certain source as sound or structure oscillations coming from dynamical wind loadings or earthquakes. Viscoelastic dampers can be designed as independent 3D parts within a mechanism but also as 2D damping treatments which consist of one or a combination of materials applied/bonded to a component to increase its ability to dissipate mechanical energy. For instance, viscoelastic dampers are useful when structures are forced to vibrate at or near their natural (resonant) frequencies, are subjected to impacts or other transient forces, or transmit vibrations to noise-radiating surfaces. The energy dissipation of a viscoelastic structure is typically quantified in terms of a loss factor, a dimensionless quantity that can be measured or predicted from a modal analysis or frequency response curves. Two categories of treatments for structural damping exist, the unconstrained layer damping (UCLD), where the material is simply attached with a strong bonding agent to the surface of a structure, and Constrained-layer damping (CLD) where the material is sandwiched between a base layer and a third constraining layer. The base layer refers to the structure to be coated. The CLD treatment provides considerably more damping effect than the free viscoelastic treatment in spite of its relative complexity. Energy dissipation is achieved in this case by shearing a viscoelastic polymer between a base structure and a constraining layer. Both types of damping structures are depicted in Figures 1 and2.

Topology optimization is essentially an iterative numerical process that seeks to find the best material layout (within a prescribed design domain) according to a given objective function and a set of design constraints, providing valuable help in problems where mechanical intuition is limited. Topology optimization of viscoelastic UCLD and CLD has been performed by many authors in the literature using different methods. Zheng et al. [START_REF] Ling | Topology optimization of constrained layer damping on plates using Method of Moving Asymptote (MMA) approach[END_REF] used the Solid Isotropic Material with Penalization (SIMP) method with the Method of Moving Asymptote (MMA) to perform topology optimization of a CLD cantilever plate treated with DYAD606 where a sum of modal loss factors is maximized. Kim [START_REF] Kim | Topology design optimization for vibration reduction: reducible design variable method[END_REF] applied SIMP and Evolutionary Structural Optimization (ESO) methods to design damping treatments for unconstrained-layer plate and shell structures. Zheng et al. [START_REF] Zheng | Topology optimization of passive constrained layer damping with partial coverage on plate[END_REF] used the same methodology to perform topology optimization of CLD with partial coverage, showing interesting performances of the optimized structure in terms of damping and mass savings. Kim et al. [START_REF] Kim | Optimal damping layout in a shell structure using topology optimization[END_REF] used the rational approximation for material properties (RAMP) with the optimality criteria method (OC) to perform topology optimization of UCLD shell structures to maximize modal loss factors. El-Sabbagh et al. [START_REF] El-Sabbagh | Topology optimization of unconstrained damping treatments for plates[END_REF] used the method presented in [START_REF] Bendsøe | Methods for optimization of structural topology, shape and material[END_REF] along with the MMA method to perform optimization of periodic and non-periodic plates. Zhanpeng et al. [START_REF] Fang | Topology Optimization for Minimizing the Resonant Response of Plates with Constrained Layer Damping Treatment[END_REF] used evolutionary structural optimization (ESO) to minimize viscoelastic CLD plate response. James et al. [START_REF] James | Topology optimization of viscoelastic structures using a time-dependent adjoint method[END_REF] used a time dependent adjoint method along with the MMA method to perform topology optimization of viscoelastically damped beams for minimum mass under time dependent loadings. Yun et al. [START_REF] Yun | Multi-material topology optimization of viscoelastically damped structures under time-dependent loading[END_REF] performed multimaterial topology optimization to maximize energy dissipation of viscoelastically damped structures subjected to unsteady loads using SIMP and MMA. Ansari et al. [START_REF] Ansari | Application of level set method to optimal vibration control of plate structures[END_REF] used a level set method to perform topology optimization of viscoelastic UCLD plate whereas van der Kolk et al. [START_REF] Van Der Kolk | Multi-material topology optimization of viscoelastically damped structures using a parametric level set method[END_REF] used a parametrized level set based method for multi-material topology optimization of beams.

The present work addresses the structural optimization of 3D fully viscoelastic structures and 2D UCLD treatments by means of the level-set method for topology optimization. For this purpose we rely on the level-set approach for multi-phase optimization detailed by Allaire et al. [START_REF] Allaire | Multi-phase structural optimization via a level set method[END_REF]. First introduced by Osher and Sethian [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF], the level-set method has the advantage of tracking the interfaces on a fixed mesh, easily managing topological changes without remeshing. Combined to the Hadamard method of shape differentiation [START_REF] Henrot | Variation et optimisation de formes: une analyse géométrique[END_REF][START_REF] Murat | Quelques résultats sur le contrôle par un domaine géométrique[END_REF][START_REF] Pironneau | Optimal shape design for elliptic systems[END_REF] in the framework of structural optimization, the level-set approach is an efficient shape and topology optimization algorithm, which gives a better description and control of the geometrical properties of the interface, avoiding typical drawbacks such as intermediate density penalization and possible spurious physical behavior during the optimization process [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF][START_REF] Wang | A level set method for structural topology optimization[END_REF]. Moreover, as remarked by Allaire et al. [START_REF] Allaire | A level-set method for vibration and multiple loads structural optimization[END_REF], the level-set method is especially well suited for vibration problems involving eigenfrequency optimization, since small holes or material islands cannot suddenly appear or disappear between two successive iterations as they do with the homogenization or SIMP methods, thus avoiding spurious modes in low density regions. Indeed, the "ersatz material approach" used within the framework of the above mentioned methods is known to produce fictitious eigenmodes localized in the weak phase which pollute the optimization process (see e.g. [START_REF] Allaire | Shape optimization by the homogenization method[END_REF], [START_REF] Bendsøe | Methods for optimization of structural topology, shape and material[END_REF], [START_REF] Bendsøe | Topology optimization: theory, methods and applications[END_REF]).

This article is organized as follows. Section 2 is devoted to the mathematical description of viscoelastic materials and the vibration problem arising for purely 3D viscoelastic structures and 2D unconstrained layer damping plates. Both cases are associated to two different models involving viscoelasticity and are dealt with in parallel inasmuch as possible. Section 3 is concerned with the study of the non-linear eigenvalue problem that stems from the aforementioned models. We provide a theoretical result characterizing the eigenvalues of such problems for a general class of viscoelastic materials. Section 4 details the sensitivity analysis of the first eigenfrequency of the structure with respect to the viscoelastic shape, and Section 5 recalls the use of the level-set method in the framework of topology optimization. Finally Section 6 outlines the numerical solution of the underlying non-linear eigenvalue problem and illustrates it with two concrete applications.

Setting of the problem

Let f be the Laplace transform of a real-valued function f (t) defined by

f (ω) = ∞ 0 e -iωt f (t)dt, ω ∈ C. (1)
From now on, we will assume a time-harmonic regime at fixed ω. For viscoelastic materials, a linear elastic constitutive relationship using Hooke's law is not an accurate model and a complex modulus is extensively used to describe the dynamic characteristics of viscoelastic materials [START_REF] Jones | Handbook of viscoelastic vibration damping[END_REF]. The stress-strain relationship of a linear viscoelastic damping material subjected to steady-state oscillatory loads thus reads:

σ(ω) = Â(ω)ˆ (ω), ( 2
)
where ω is the complex pulsation, Â is the complex (or dynamic) elastic tensor and σ and ˆ are the Laplace transforms of stress and the strain, respectively. This relation stems from the Laplace transform of the stress history given by

σ(t) = t 0 Y (t -τ ) d (τ ) dτ dτ, (3)
with Y the relaxation function (or relaxation modulus) and

 = iω Ŷ (ω). ( 4
)
The relaxation function Y accounts for the material stress behavior in a relaxation test (gradual disappearance of stresses from a viscoelastic medium after it has been deformed) and the Boltzmann superposition principle (the state of stress or deformation of a viscoelastic body is a function of all the stresses and/or deformations applied to the material) [START_REF] Findley | Creep and relaxation of nonlinear viscoelastic materials[END_REF]. In other words, the relaxation function Y governs the decrease of the stress towards an asymptotical limit in time when the viscoelastic material is subjected to an instantaneous constant strain at t = 0. Then for a sequence d (t) of applied strains in time, the stress is simply equal to the "sum of products" of the delayed relaxation function and the constant applied strains. This integral is later expressed rather as a convolution.

For the sake of simplicity, we will henceforth write (A, σ, ) instead of ( Â, σ, ˆ ) unless we want to explicitly emphasize the Laplace transform.

Let (Ω, O) be two bounded open sets with Ω ⊂ O ⊂ R d (d ∈ {2, 3}).

• If d = 3, we will say that Ω is occupied by a viscoelastic material with complex elastic tensor A and density ρ > 0 within a reference (or working) domain O. We will call this configuration the 3D viscoelastic structure. The boundary of Ω is made of two disjoint parts

∂Ω = Γ N ∪ Γ D ,
with Dirichlet boundary conditions on Γ D and Neumann boundary conditions on Γ N .

• If d = 2, O will be the reference configuration for a plate in the xy-plane and Ω the surface occupied by the viscoelastic coating. We will call this configuration the composite sandwich structure. By abuse of notation, we split the boundary of O (∂O = Γ N ∪ Γ D ) instead of Ω since no boundary conditions are imposed on ∂Ω.

3D viscoelastic structure

We denote by ω ∈ C the complex pulsation and by û ∈ H 1 (Ω; C) 3 the associated mode, i.e. the corresponding displacement field in Ω. For the sake of simplicity, we will generally write u unless we want to explicitly emphasize the Laplace transform. The pair (ω, u) is a solution of the non-linear eigenvalue problem

   -div(A(ω)e(u)) = ω 2 ρ u in Ω, u = 0 on Γ D , A(ω)e(u) • n = 0 on Γ N , (5)
where e(u) = ∇u + ∇u T 2 is the symmetrized strain tensor and n is the normal vector. The state equation in [START_REF] El-Sabbagh | Topology optimization of unconstrained damping treatments for plates[END_REF] stems from applying the constitutive equation ( 2) and the Laplace transform to the evolution (wave) equation of linear elasticity

ρu tt -div(σ) = 0 in Ω. ( 6
)
We remark that the eigenvalue problem ( 5) is non-linear since the complex elastic (or stiffness) tensor A also depends on ω.

Composite sandwich structure

For modeling vibration damping treatments we use a plate model approximation of [START_REF] El-Sabbagh | Topology optimization of unconstrained damping treatments for plates[END_REF]. Consider the small transverse vibration of a uniform plate with thickness h(x), density ρ(x) and reference configuration for the plate O with d = 2. The transverse displacement of x at time t is denoted as w(x, t) but for the sake of simplicity, we rather denote as w (instead of ŵ) the vertical displacement mode associated with the complex pulsation ω ∈ C. Unconstrained-layer damping (UCLD) is one of the simplest forms of material application. The material is simply attached with a strong bonding agent to the surface of a structure. Energy is dissipated as a result of extension and compression of the damping material under flexural stress from the base structure. From now on, for every general quantity ξ defined either on the viscoelastic coating, the base plate or the global composite (viscoelastic coating plus base plate), we will use the notations ξ c , ξ p , ξ g respectively. The following assumptions are made:

• Transverse shear and rotational and in-plane inertia effects in both the plate and the coating are negligible for the lower modes of vibration,

• There are no applied in-plane loads,

• Displacements are small and changes in thickness are negligible,

• The viscoelastic coating is applied to only one side of the plate

• The plate and the coating are homogeneous and isotropic, and subjected to a state of plane stress,

• Displacements are continuous across the interface between the plate and the coating,

• The Poisson's ratio of the coating is a real constant and the coating is incompressible (ν = 0.5).

Under these assumptions, classical Kirchhoff-Love plate theory can be used [START_REF] Reddy | Mechanics of laminated composite plates: theory and analysis[END_REF]. It is worth mentioning that the plate is considered to be in a free vibration state. Let the pair (ω, w), with w ∈ H 2 (O; C), be the solution of the non-linear eigenvalue problem

   -∇ 2 : (D(ω)∇ 2 w) = ω 2 ρ w in O, w = 0, ∇w = 0 on Γ D , n • (D(ω)∇ 2 w) • n = 0 on Γ N . (7)
Neutral plane of base plane

Neutral plane of the viscoelastic layer

Neutral plane of the composite where ∇ 2 is the second order partial derivatives operator. The aforementioned state equation follows from the virtual work principle applied to the total energy of the plate (bending strain energy plus kinetic energy) and the Laplace transform. As mentioned before, let Ω be the region of O where the coating is applied.

h c h p h n h pc z = 0 z x y
Then the function ρ(x) is piecewise constant and corresponds to the surface density

ρ(x) = ρ p := ρ V p h p if x ∈ O\Ω, ρ g := (ρ V p h p + ρ V c h c ) if x ∈ Ω, with ρ V p , ρ V c
and h p , h c the volumetric densities and thicknesses of the plate and the viscoelastic coating, respectively. The equivalent bending stiffness tensor D for an isotropic base and coating material reads [START_REF] Parthasarathy | Partial coverage of rectangular plates by unconstrained layer damping treatments[END_REF]

D(ω, x) =        D p := h 3 p 12 A p if x ∈ O\Ω, D g := h 3 p 12 + h p h 2 n A p + h 3 c 12 + h c (h pc -h n ) 2 A c (ω) if x ∈ Ω, (8)
where the tensors A p , A c (ω) stand for the extensional stiffness tensors of the base and viscoelastic materials respectively (we remark that only the coating material depends on ω, more specifically its Young modulus),

A l = E l (1 -ν 2 l )   1 ν l 0 ν l 1 0 0 0 0.5(1 -ν l )   , with the index l = p, c, (9)
h pc represents the distance between the neutral planes of the base plate and the viscoelastic layer material and h n (see Figure 3) represents the distance between the composite neutral axis and the base plate neutral plane which can be approximated as

h n = E c h c h pc E p h p + E c h c . ( 10
)
The structure of the tensor D ensues from the integration of the laminate through the thickness moment resultants meanwhile the correction terms h n , h pc arise when h p = h c in order to compensate for the shift of the neutral fiber of the laminate, i.e. the place in the plate where the thickness-averaged stresses vanish [START_REF] Parthasarathy | Partial coverage of rectangular plates by unconstrained layer damping treatments[END_REF].

Optimization problem

For the sake of clarity, we detail the optimization problem only within the damping optimization context exposed in Section 2.1. The same notions however remain valid for viscoelastic damping treatments. Supposing that (5) admits a countably infinite family of solutions (ω k , u k ) k≥1 in C × H 1 (Ω; C) 3 (see more details in the next section), with the eigenfunctions or modes normalized by imposing that Ω ρ|u k | 2 dx = 1, we are interested in solving the optimization problem sup

Ω∈U ad η(Ω), (11)
where

U ad = {Ω ⊂ O such that |Ω| = V} (12)
is the set of admissible shapes with fixed volume V, O is the working domain (a bounded subset of R 3 ) and η(Ω) is the modal loss factor of the structure for its first eigenvalue:

η(Ω) := Im(ω 2 1 ) Re(ω 2 1 ) = Im Ω Ae(u 1 ) : e(u 1 )dx Re Ω Ae(u 1 ) : e(u 1 )dx , ( 13
)
where u is the complex conjugate of u. The value of η represents the ratio between dissipated and stored energy or in other words, the structure capacity to dissipate energy. It is obviously a measure of the structure mechanical damping [START_REF] Carfagni | The loss factor as a measure of mechanical damping[END_REF].

Remark 2.1. It is well-known that ( 11) is usually not well posed on the set of admissible shapes (12) (i.e. it has no solution). In order to obtain existence of optimal shapes, some smoothness or geometrical or topological constraints are required (see for instance [START_REF] Ambrosio | An optimal design problem with perimeter penalization[END_REF][START_REF] Chambolle | A density result in two-dimensional linearized elasticity, and applications[END_REF][START_REF] Chenais | On the existence of a solution in a domain identification problem[END_REF] ). In the present article, we shall work with a subset of ( 12) composed of smoother shapes in order to define properly the notion of shape derivative (consult Section 4).

Similarly to the notion of material loss factor of a viscoelastic material, which allows to express the dissipative behavior of the complex Young and shear moduli E(ω), G(ω) via the Young and shear moduli factors η E , η G :

E = Re(E)(1 + iη E ) and G = Re(G)(1 + iη G ), (14)
the concept of modal loss factor accounts for the dissipative behavior of the structure first resonant frequency

ω 2 1 = Re(ω 2 1 )(1 + iη).
Since the eigenvalues of ( 5) cannot be naturally ordered in C, ω 1 will be defined as the closest eigenvalue of (5) to ω 1 (in the sens of modulus 3 and ω k ≤ ω k+1 for all k ≥ 1, is the countably infinite family of solutions of the purely elastic problem associated with (5):

|ω 1 -ω 1 |), where (ω k , u k ) k≥1 in R × H 1 (Ω; R)
-div(A(0)e(u)) = ω 2 ρ u in Ω. ( 15
)

The non-linear eigenvalue problem

Problem (5) can be cast as a generalized eigenvalue problem

T (ω)u = 0 ( 16
)
where T (ω) is a linear operator depending (non-linearly) on a parameter ω. In particular for ( 5)

T (ω) = -div(A(ω)e(•)) -ω 2 ρ Id(•),
where Id() is the identity operator. A solution u = 0 will exist only for some particular values of ω (also called eigenvalues).

We devote this section to the study of the characterization of the solutions of ( 5) and [START_REF] Fang | Topology Optimization for Minimizing the Resonant Response of Plates with Constrained Layer Damping Treatment[END_REF]. We prove in particular that the underlying spectrum of the aforementioned problems is discrete. This feature is of paramount importance since a continuous spectrum or an accumulation point in the spectrum may lead to spurious eigenvalues in numerical applications.

Preliminary results

Let V and W be two complex Hilbert spaces with V compactly embedded in W. Let also {E ω } ω∈D , with D ⊂ C open bounded, be a family of complex mappings E ω : V × V → C. Definition 3.1 (Sesquilinear form). We say that {E ω } ω∈D is a family of sesquilinear forms if for each ω

1. E ω (u 1 + u 2 , v 1 + v 2 ) = E ω (u 1 , v 1 ) + E ω (u 1 , v 2 ) + E ω (u 2 , v 1 ) + E ω (u 2 , v 2 ), 2. E ω (αv 1 , βv 2 ) = αβE ω (v 1 , v 2 ), for all u 1 , u 2 , v 1 , v 2 ∈ V and all α, β ∈ C.
Definition 3.2 (Continuous form). We say that {E ω } ω∈D is a family of continuous forms if there exists a function d :

D → R + such that |E ω (v 1 , v 2 )| ≤ d(ω) v 1 V v 2 V , ∀v 1 , v 2 ∈ V. ( 17
)
We will suppose additionally that the function d(ω) is continuous in D (thus bounded).

Definition 3.3 (Coercive form). We say that {E ω } ω∈D is a family of coercive forms if there exists a function c :

D → R + such that |E ω (v, v)| ≥ c(ω) v 2 V , ∀v ∈ V. ( 18
)
We will suppose additionally that the function c(ω) is continuous and c(ω) ≥ α > 0 in D (thus bounded away from 0). Definition 3.4 (Holomorphic form). We say that {E ω } ω∈D is a family of holomorphic forms if for every fixed pair

v 1 , v 2 ∈ V the mapping ω ∈ C → E ω (v 1 , v 2 ) ∈ C
is differentiable at ω for every point ω in D.

From now on, we will assume that {E ω } ω∈D is a family of continuous coercive sesquilinear forms with continuous bounds.

Proposition 3.1. For each ω ∈ D, define the operators G ω , G † ω : W → W such that for given g ∈ W, G ω (g) and G † ω (g) are respectively the solutions of

E ω (G ω (g), v) = g, v W , ∀v ∈ V, (19a)
E ω (v, G † ω (g)) = g, v W , ∀v ∈ V. (19b)
Then for each ω ∈ D, the operators G ω , G † ω are well-defined and compact.

Proof. We only give the proof for G ω as the argument is analogous for G † ω . According to ( 17) and ( 18), we can apply the Lax-Milgram theorem to (19a) so for a given g ∈ W, there exists a unique solution

G ω (g) ∈ V ⊂ W. Furthermore, since V is compactly embedded in W, there exists a constant b > 0 such that b G ω (g) W ≤ G ω (g) V , so the operator G ω is continuous with G ω (g) W ≤ 1 b 2 α g W , (20)
thanks to Cauchy-Schwarz inequality and the Lax-Milgram estimate of (19a).

Proposition 3.2. Assume that {E ω } ω∈D is a family of holomorphic forms and for ω fixed, the first and second derivative mappings of E ω , denoted respectively as

∂ ω E ω : V × V → C, ∂ 2 ω E ω : V × V → C,
are continuous sesquilinear forms with continuous bounds. For every ω ∈ D, define the complex mapping

F ω : W × W → C as F ω (v 1 , v 2 ) = G ω (v 1 ), v 2 W , v 1 , v 2 ∈ W, (21)
where the operator G ω is defined in (19a). Then {F ω } ω∈D is a family of holomorphic forms.

Proof. We will prove that for an arbitrary fixed pair v 1 , v 2 ∈ W, the mapping

ω ∈ D → F ω (v 1 , v 2 ) = G ω (v 1 ), v 2 W ∈ C (22)
is differentiable at ω 0 for every ω 0 ∈ D. Let p ω0 ∈ V be the solution of

E ω0 (p ω0 , v) = -∂ ω E ω0 (G ω0 (v 1 ), v), ∀v ∈ V ( 23
)
for a given v 1 ∈ V. The existence and uniqueness of p ω0 stem from the Lax-Milgram theorem according to [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF], [START_REF] Wang | A level set method for structural topology optimization[END_REF] and the continuity of the operator ∂ ω E ω0 . Let ∆ω ∈ C be small enough so ω 0 + ∆ω ∈ D. Then ( 22) is differentiable at

ω 0 with lim ∆ω→0 G ω0+∆ω (v 1 ), v 2 W -G ω0 (v 1 ), v 2 W ∆ω -p ω0 (v 1 ), v 2 W = 0. (24)
Indeed

G ω0+∆ω (v 1 ), v 2 W -G ω0 (v 1 ), v 2 W ∆ω -p ω0 (v 1 ), v 2 W = 1 |∆ω| G ω0+∆ω (v 1 ) -G ω0 (v 1 ) -∆ω p ω0 (v 1 ), v 2 W = 1 |∆ω| v 2 , G ω0+∆ω (v 1 ) -G ω0 (v 1 ) -∆ω p ω0 (v 1 ) W = 1 |∆ω| E ω0+∆ω G ω0+∆ω (v 1 ) first member -G ω0 (v 1 ) -∆ω p ω0 (v 1 ) second member , G † ω0+∆ω (v 2 ) . (25)
Then developing each member of ( 25) :

E ω0+∆ω G ω0+∆ω (v 1 ), G † ω0+∆ω (v 2 ) = v 1 , G † ω0+∆ω (v 2 ) W (26)
according to (19a) and 28)

-E ω0+∆ω G ω0 (v 1 ) + ∆ω p ω0 (v 1 ), G † ω0+∆ω (v 2 ) = -E ω0 G ω0 (v 1 ) + ∆ω p ω0 (v 1 ), G † ω0+∆ω (v 2 ) ( 27
)
-∆ω∂ ω E ω0 G ω0 (v 1 ) + ∆ω p ω0 (v 1 ), G † ω0+∆ω (v 2 ) (
-(∆ω) 2 ∂ 2 ω E ω0 G ω0 (v 1 ) + ∆ω p ω0 (v 1 ), G † ω0+∆ω (v 2 ) (29) + O((∆ω) 3 v 1 V v 2 V )
according to the second-order Taylor expansion of E ω with respect to ω (by hypothesis we assume the mapping ω → E ω (v 1 , v 2 ) holomorphic for ω ∈ D). By virtue of (19a) and ( 23), [START_REF] Ambrosio | An optimal design problem with perimeter penalization[END_REF] satisfies

-E ω0 G ω0 (v 1 ) + ∆ω p ω0 (v 1 ), G † ω0+∆ω (v 2 ) = -v 1 , G † ω0+∆ω (v 2 ) W + ∆ω∂ ω E ω0 (G ω0 (v 1 ), G † ω0+∆ω (v 2 )),
so adding ( 26), ( 27), ( 28) and ( 29), a few terms cancel and we deduce that ( 25) is equal to

= 1 |∆ω| -(∆ω) 2 ∂ ω E ω0 p ω0 (v 1 ), G † ω0+∆ω (v 2 ) + ∂ 2 ω E ω0 G ω0 (v 1 ), G † ω0+∆ω (v 2 ) + O((∆ω) 3 v 1 V v 2 V ) ≤ O(|∆ω| v 1 V v 2 V ),
where the last inequality stems from boundedness in D of the continuity bounds of ∂ ω E ω , ∂ 2 ω E ω and the Lax-Milgram bounds derived in the proof of Proposition 3.1 (reasoning by analogy, the same bounds can be applied to ( 23)). Finally taking the limit when ∆ω → 0 yields (24).

Main result

The main conclusion of this section corresponds to the application of a result of standard analytical perturbation theory of compact operators. This theorem describes locally the eigenvalues of a general non-linear eigenvalue problem Theorem 3.1 ( [START_REF] Kato | Perturbation theory for linear operators[END_REF], chapter VII, Th. 1.9). Let S(ω) be a family of compact operators (for each ω fixed) and holomorphic with respect to ω ∈ D ⊂ C bounded. Define ω as a singular point if 1 is an eigenvalue of S(ω). Then either all ω ∈ D are singular points or there are only a finite number of singular points in each compact subset of D. Theorem 3.2. Suppose that the weak formulation of (16) has the following structure on ω

E ω (u, v) -ω 2 u, v W = 0, ∀v ∈ V, (30)
and assume the same hypotheses as Proposition 3.1 and Proposition 3.2. Furthermore suppose that 0 ∈ D and E 0 hermitian, i.e.

E 0 (u, v) = E 0 (v, u), ∀u, v ∈ V.
Then there are only a finite number of eigenvalues of (30) in D (i.e. values of ω for which there exist u = 0 solution of (30)).

Proof. Define the operator S ω : W → W as

S ω = ω 2 G ω
Then according to Proposition 3.1 and Proposition 3.2, {S ω } ω∈C is a family of compact operators and S ω is holomorphic with respect to ω in D. Furthermore, ω = 0 is an eigenvalue of (30) if and only if 1 is an eigenvalue of S ω according to (19a). If every ω ∈ D\{0} fulfilled [START_REF] Kato | Perturbation theory for linear operators[END_REF], then since D is open, there would be a sequence {z j } j∈N ∈ D and {v j } j∈N ∈ V with v j W = 1, such that z j → 0 when j → ∞ and

|E zj (v j , v j )| = |z 2 j v j , v j W | ≤ |z 2 j |, ∀j ∈ N. (31)
Now since E 0 is hermitian and coercive (the family of operators {E ω } ω∈D was supposed coercive from the beginning), thus for all v = 0, there exists λ > 0 such that

E 0 (v, v) v V ≥ λ, so 0 < λ = λ v j 2 W ≤ λ v j 2 V ≤ E 0 (v j , v j ), ∀j, (32)
since V is embedded in W. But then, using the fact that E ω is holomorphic and (31)

E 0 (v j , v j ) ≤ |E 0 (v j , v j ) -E zj (v j , v j )| + |E zj (v j , v j )| = O(|z j |)
and we obtain a contradiction with [START_REF] Murat | Etude de problèmes d'optimal design[END_REF] for j large enough. The desired result finally stems from Theorem 3.1.

We give an application of the above theorem to the 3D viscoelastic and UCLD plate vibration problems Proof. The result derives from Theorem 3.2 by considering

E ω (u, v) = Ω A(ω)e(u) : e(v)dx, V = v ∈ H 1 (Ω; C) 3 : v = 0 on Γ D , W = L 2 (Ω; C) 3
in problem ( 5) and

E ω (w, ζ) = O D(ω)∇ 2 w : ∇ 2 ζdx, V = ζ ∈ H 2 (O; C) : ζ = 0, ∇ζ = 0 on Γ D , W = L 2 (O; C)
in problem [START_REF] Fang | Topology Optimization for Minimizing the Resonant Response of Plates with Constrained Layer Damping Treatment[END_REF].

4. Sensitivity analysis of ω 1 with respect to the shape Ω

In order to apply a gradient method to the minimization of (11), we dedicate this section to the study of the shape differentiability of ω 1 and eigenvectors (u 1 , w 1 ), solutions of ( 5) and [START_REF] Fang | Topology Optimization for Minimizing the Resonant Response of Plates with Constrained Layer Damping Treatment[END_REF] respectively. Nevertheless all results remain valid for k ≥ 2 by replacing (ω 1 , u 1 , w 1 ) with (ω k , u k , w k ).

Before going any further, we recall the classical notion of shape derivative which goes back, at least, to Hadamard (see the modern reference books [START_REF] Pironneau | Optimal shape design for elliptic systems[END_REF], [START_REF] Soko Lowski | Introduction to shape optimization: Shape sensitivity analysis[END_REF]). Here, we follow the approach of Murat and Simon [START_REF] Murat | Etude de problèmes d'optimal design[END_REF], [START_REF] Simon | Differentiation with respect to the domain in boundary value problems[END_REF] . Starting from a smooth reference open set Ω , we consider domains of the type

Ω θ = (Id + θ)(Ω),
where Id is the identity mapping from R d into R d , and θ ∈ W 1,∞ (R d ; R d ). It is well-known that, for sufficiently small θ, (Id + θ) is a diffeomorphism in R d . We denote as C(Ω) the set of admissible shapes defined as

{(Id + θ)(Ω), θ ∈ W 1,∞ (R d ; R d )}.
Definition 4.1 (Shape derivative). Let J(Ω) be a mapping from the set of admissible shapes C(Ω) to R. The shape derivative of J(Ω) at Ω is defined as the Fréchet derivative in W 1,∞ (R d ; R d ) at 0 of the mapping θ → J((Id + θ)(Ω)), i.e.,

J((Id

+ θ)(Ω)) = J(Ω) + J (Ω)(θ) + o(θ) with lim θ→0 |o(θ)| |θ| = 0, ( 33
)
where

J (Ω) is a continuous linear form on W 1,∞ (R d ; R d ).
In his research on elastic plates [START_REF] Hadamard | Mémoire sur le problème d'analyse relatif à l'équilibre des plaques élastiques encastrées[END_REF] in 1907, Hadamard showed how to obtain the shape derivative of a functional J(Ω) by considering normal perturbations of the boundary Ω. This fundamental result was made rigorous later by Zolésio [START_REF] Delfour | of advances in design and control[END_REF] in the so-called "structure theorem" which states that when J(Ω) and Ω are smooth enough, one may write the shape derivative of J(Ω) as an integral over ∂Ω, i.e. there exists V ∈ L 1 (∂Ω) such that

J (Ω)(θ) = ∂Ω (θ • n)Vds. (34)
To fix ideas, we give two examples of shape derivative.

(x) ∈ W 1,1 (R d ). Define J 1 (Ω) = Ω φ(x)dx.
Then J 1 is differentiable at Ω and

J 1 (Ω)(θ) = ∂Ω θ(x) • n(x)φ(x)ds for any θ ∈ W 1,∞ (R d ; R d ). For φ(x) ∈ W 2,1 (R d ) define J 2 (Ω) = ∂Ω φ(x)ds.
Then J 2 is differentiable at Ω and

J 2 (Ω)(θ) = ∂Ω θ(x) • n(x) ∂φ ∂n + Hφ ds, for any θ ∈ W 1,∞ (R d ; R d ),
where H is the mean curvature of ∂Ω defined by H = div n Before computing the shape derivatives of our damping problem, we introduce the definitions of the Lagrangian and Eulerian derivatives. Definition 4.2. The function u : Ω → u(Ω) ∈ H 1 (Ω) admits a Lagrangian (or material) derivative Y (Ω)(θ) at a given domain Ω, provided that the application

θ ∈ W 1,∞ (R d ; R d ) → ũθ (Ω) := u((Id + θ)Ω) • (Id + θ) ∈ H 1 (Ω).
is Fréchet differentiable at θ = 0. Then Y (Ω 0 )(θ) = ũ θ (Ω). Definition 4.3. The function u : Ω → u(Ω) admits a Eulerian derivative u (Ω)(θ) at a given domain Ω, if u admits a Lagrangian derivative Y (Ω)(θ) at Ω in the direction θ and ∇u(Ω) • θ ∈ H 1 (Ω). In that case

u (Ω)(θ) = Y (Ω)(θ) -∇u(Ω) • θ.

Differentiability of ω 1

Theorem 4.1. Suppose Ω measurable and bounded and ω 1 (Ω) simple. Then the first eigenfrequency ω 1 (Ω), u 1 (Ω) (for d = 3) and w 1 (Ω) (for d = 2) are shape differentiable.

Proof. We give a sketch of proof since it corresponds to an adaptation of the result exhibited in [START_REF] Henrot | Variation et optimisation de formes: une analyse géométrique[END_REF] chapter 5, page 210. Also, as a matter of simplicity, we detail the proof argument only for d = 3.

Define u θ and ω(θ) as the solutions of the problem [START_REF] El-Sabbagh | Topology optimization of unconstrained damping treatments for plates[END_REF] within Ω θ = (Id + θ)Ω, where θ ∈ W 1,∞ (R 3 ; R 3 ). In other words, -div(A(ω)e(u θ ))ω 2 ρu θ = 0 in Ω θ ,

Ω θ ρ|u θ | 2 dx -1 = 0. ( 35
)
Let us consider the operator [START_REF] De Gournay | Velocity extension for the level-set method and multiple eigenvalues in shape optimization[END_REF], such that F(θ, u, ω) = 0 represents the transported problem [START_REF] Delfour | of advances in design and control[END_REF] on Ω. The main idea of the proof is to apply the implicit function theorem to F in order to show that there exist differentiable functions ω 1 (θ) and u 1 (θ) at θ = 0, solutions of [START_REF] Delfour | of advances in design and control[END_REF]. For that purpose, we need to prove that the differential operator D v,ω F (0, u 1 , ω 1 ) given by D v,ω F(0, u 1 , ω 1 )(ṽ, μ) =div(A(ω 1 )e(ṽ))ω 2 1 ρṽμ(div(∂ ω A(ω 1 )e(u 1 )) + 2ω 1 ρu 1 ), 2

F : W 1,∞ (R 3 ; R 3 ) × V × C → H -1 (Ω; C 3 ) × R, with V defined by
Ω ρṽ • ū1 dx ,
for every (ṽ, μ) ∈ V × C is an isomorphism (by virtue of the inverse function theorem in Banach spaces) from

V × C on H -1 (Ω; C 3 ) × C. Since D v,µ F(0, u 1 , ω 1
) is continuous, according to the Banach theorem, to prove that D v,µ F(0, u 1 , ω 1 ) is an isomorphism one simply needs to show that D v,µ F(0, u 1 , ω 1 ) is a bijection. Thus, given (Z, Λ) ∈ H -1 (Ω; C 3 ) × C, we need to establish that the problem -div(A(ω 1 )e(ṽ))ω 2 1 ρṽμ(div(∂ ω A(ω 1 )e(u 1 )) + 2ω

1 ρu 1 ) = Z in Ω, 2 Ω ρṽ • ū1 dx = Λ
has an unique solution (ṽ, μ) ∈ V × C. This is true thanks to the fact that for a fixed ω ∈ C (in our case ω = ω 1 ), the same argument as that given in [START_REF] Henrot | Variation et optimisation de formes: une analyse géométrique[END_REF] can be applied, namely in a variational sens the operator

(-div(A(ω)e(•))) -1 : V → V ⊂ V ,
with V the dual of V , is compact and thus the Fredholm alternative theorem can be applied to the operator

-div(A(ω 1 )e(•)) -ω 2 1 ρ Id(•)
whose kernel is of dimension one since we assume ω 1 simple. Thus, in order to be unique the solution of the above system, the following orthogonality condition must stand Zμ(div(∂ ω A(ω 1 )e(u 1 )) + 2ω 1 ρu 1 ), u 1 = 0, which determines uniquely the value of μ. The uniqueness of ṽ stems from the decomposition ṽ = v 0 + su 1 , where v 0 is a particular solution and s ∈ R. s is uniquely determined by the second condition of the system.

Shape derivative in the case of a 3D viscoelastic structure

Theorem 4.2. Define the space

V = v ∈ H 1 (Ω; C) 3 : v = 0 on Γ D . (36)
Let (ω 1 , u 1 ) ∈ C × V be the solution of the eigenvalue problem (5) and ω 1 simple. Also let A H (ω 1 ) = A(ω 1 ) T be the conjugate transpose tensor of A(ω 1 ) and introduce the adjoint eigenvector p 1 ∈ V solution of

Ω A H (ω 1 )e(p 1 ) : e(v)dx = ω1 2 Ω ρ(p 1 • v)dx, ∀v ∈ V, (37)
with Ω ρ|p 1 | 2 dx = 1 and ω2 1 = ω 2 1 . Then:

ω 1 (Ω)(θ) = -∂Ω θ • n ω 2 1 ρ (u 1 • p 1 ) -A(ω 1 )e(u 1 ) : e(p 1 ) ds Ω 2ω 1 ρ (u 1 • p 1 ) -∂ ω A(ω 1 )e(u 1 ) : e(p 1 ) dx . ( 38
)
Proof. As mentioned in the previous section, the fact that ω 1 is simple ensures that ω 1 (Ω) and u 1 (Ω) are shape-differentiable. Thus we can differentiate with respect to the domain both sides of the variational formulation of ( 5)

Ω ρ(u 1 • v)dx, ∀v ∈ V (39) Taking θ ∈ W 1,∞ (R 3 , R 3 ) with θ = 0 on Γ N ∪ Γ D : ∂Ω A(ω 1 )e(u 1 ) : e(v)(θ • n)ds + ω 1 (Ω)(θ) Ω (∂ ω A)e(u 1 ) : e(v)dx + Ω A(ω 1 )e(u 1 (Ω)(θ)) : e(v)dx = 2ω 1 ω 1 (Ω)(θ) Ω ρ(u 1 • v)dx + ω 2 1 Ω ρ(u 1 (Ω)(θ) • v)dx + ω 2 1 ∂Ω ρ(u 1 • v)(θ • n)ds (40)
Applying the conjugate on both sides of equation ( 37) one gets

Ω A(ω 1 )e(v) : e(p 1 )dx = ω 2 1 Ω ρ(v • p1 )dx, ∀v ∈ V, (41)
so taking v = p 1 in ( 40) and v = u 1 (Ω)(θ) in ( 41), one realizes that all the terms containing u 1 (Ω)(θ) in ( 40) cancel. Thus [START_REF] Johnson | Finite element prediction of damping in structures with constrained viscoelastic layers[END_REF] yields

∂Ω A(ω 1 )e(u 1 ) : e(p 1 )(θ • n)ds + ω 1 (Ω)(θ) Ω (∂ ω A)e(u 1 ) : e(p 1 )dx = 2ω 1 ω 1 (Ω)(θ) Ω ρ(u 1 • p 1 )dx + ω 2 1 ∂Ω ρ(u 1 • p 1 )(θ • n)ds,
from which (38) follows.

Remark 4.1. When the tensor A is independent from the frequency ω, we recover the well-known formula [START_REF] Allaire | A level-set method for vibration and multiple loads structural optimization[END_REF] for the derivative:

ω 1 (Ω)(θ) = -∂Ω θ • n ω 2 1 ρ (u 1 • p 1 ) -A(ω 1 )e(u 1 ) : e(p 1 ) ds Ω 2ω 1 ρ (u 1 • p 1 )dx
Furthermore if A is real (since by definition is symmetric) then the problem becomes self-adjoint (i.e. p 1 = u 1 ).

Shape derivative of a composite sandwich structure

We refer to the notations given in Section 2.2. Define a bounded working domain O ⊂ R 2 and an admissible shape Ω ⊂ O occupied by a viscoelastic damping treatment. Within this framework, ( 7) can be formulated as a multi-phase problem. As mentioned in [START_REF] Allaire | Multi-phase structural optimization via a level set method[END_REF], a sharp interface between the different phases leads to numerical difficulties that can be avoided by considering a smeared or diffuse interface approach. We follow herein this smoothed formulation and as abuse of notation we consider the same eigenelements ω 1 and w 1 . Let D and ρ be the smooth interpolations between D g , D p , (with respective surface densities ρ g , ρ p ) occupying Ω and O\Ω respectively, defined by

D = D g -H ε (d Ω )∆D, ∆D := (D g -D p ), ρ = ρ g -H ε (d Ω )∆ρ, ∆ρ := (ρ g -ρ p ), ( 42
)
where H ε (r) : R → [0, 1] is a smooth approximation of the Heaviside function

H ε (r) =    0 if r < -ε, 1 2 (1 + r ε + 1 π sin πr ε ) if -ε ≤ r ≤ ε, 1 if r > ε, (43)
and d Ω (x) : R d → R the signed distance function at x associated with Ω

d Ω (x) =      0 if x ∈ ∂Ω, -min x I ∈∂Ω |x -x I | if x ∈ Ω, min x I ∈∂Ω |x -x I | if x ∈ O\Ω.
Before giving the shape derivative of ω 1 in the context of Section 2.2, we collect some general definitions and results related to the signed distance function d Ω in R d (for more details consult [START_REF] Allaire | Multi-phase structural optimization via a level set method[END_REF]). • For any x ∈ R d we define the set of projections of x on ∂Ω as

Π ∂Ω (x) = {y 0 ∈ ∂Ω such that |x -y 0 | = inf y∈∂Ω |x -y|}.
When Π ∂Ω (x) reduces to a single point, we will call it the projection P ∂Ω (x) of x onto ∂Ω.

• The skeleton Σ of ∂Ω is defined as

Σ := {x ∈ R d such that (d Ω ) 2 is not differentiable at x}. Lemma 4.2. Let Ω ⊂ R d be a Lipschitz bounded open set. • A point x /
∈ ∂Ω has a unique projection P ∂Ω (x) on ∂Ω if and only if x / ∈ Σ.

• Σ has zero Lebesgue measure in R d .

The signed distance function is not, strictly speaking, shape differentiable in the sense of [START_REF] Simon | Differentiation with respect to the domain in boundary value problems[END_REF]. One reason is the lack of smoothness of the gradient of d Ω at the skeleton Σ. However, its pointwise values d Ω (x) are shape differentiable for x ∈ O\Σ. In particular for d = 2 (the result also holds for d = 3) Proposition 4.1. Assume Ω ⊂ O is an open set of class C 1 , and fix a point x / ∈ Σ. Then θ → d (Id+θ)Ω (x) is Gâteaux-differentiable at θ = 0, as a mapping from W 1,∞ (O; R 2 ) into R, and its derivative is

d Ω (θ)(x) = -θ(P ∂Ω (x)) • n(P ∂Ω (x)). ( 44
) Theorem 4.3. Define the space W = ζ ∈ H 2 (O; C) : ζ = 0, ∇ζ = 0 on Γ D . Let (ω 1 , w 1 )
∈ C × W be the solutions of the eigenvalue problem (7) and assume ω 1 simple. Also let D H (ω 1 ) = D(ω 1 ) T be the conjugate transpose tensor of D and introduce the adjoint eigenvector

q 1 ∈ W solution of O D H (ω 1 )∇ 2 q : ∇ 2 ζdx = ω2 O ρ(q • ζ)dx, ∀ζ ∈ W, (45)
with O ρ|q 1 | 2 dx = 1 and ω2 1 = ω 2 1 . Then:

ω 1 (Ω)(θ) = -O ∂ r H (d Ω ) ω 2 1 ∆ρ(w 1 • q 1 ) -∆D(ω 1 )∇ 2 w 1 : ∇ 2 q 1 θ(P ∂Ω (x)) • n(P ∂Ω (x))dx O 2ω 1 ρ (w 1 • q 1 ) -∂ ω D(ω 1 )∇ 2 w 1 : ∇ 2 q 1 dx , (46)
where P ∂Ω (x) is the projection of x onto ∂Ω.

Proof. The shape differentiability of ω 1 (Ω) follows from the same argument as in the proof of Theorem 4.1.

Hence we differentiate with respect to the domain both sides of the variational formulation of (5) but this time within the fixed domain

O O D(ω 1 )∇ 2 w 1 : ∇ 2 ζdx = ω 2 1 O ρ(w 1 • ζ)dx, ∀ζ ∈ W. Taking θ ∈ W 1,∞ (R 3 , R 3 ) with θ = 0 on ∂O O D (Ω)(θ)(ω 1 )∇ 2 w 1 : ∇ 2 ζdx + ω 1 (Ω)(θ) O ∂ ω D∇ 2 w 1 : ∇ 2 ζdx + O D(ω 1 )∇ 2 w 1 (Ω)(θ) : ∇ 2 ζdx = 2ω 1 ω 1 (Ω)(θ) O ρ(w 1 • ζ)dx + ω 2 1 O ρ(w 1 (Ω)(θ) • ζ)dx + ω 2 1 O ρ (Ω)(θ)(w 1 • ζ)dx. ( 47
)
Applying the same argument given in the proof of Theorem 4.2, by choosing the right test functions in (47) (ζ = q 1 ) and ( 45)

(ζ = w 1 (Ω)(θ)), (47) yields
O D (Ω)(θ)(ω 1 )∇ 2 w 1 : ∇ 2 q 1 dx + ω 1 (Ω)(θ) O ∂ ω D∇ 2 w 1 : ∇ 2 q 1 dx = 2ω 1 ω 1 (Ω)(θ) O ρ(w 1 • q 1 )dx + ω 2 1 O ρ (Ω)(θ)(w 1 • q 1 )dx, so ω 1 (Ω)(θ) = -O ω 2 1 ρ (Ω)(θ)(w 1 • q 1 ) -D (Ω)(θ)(ω 1 )∇ 2 w 1 : ∇ 2 q 1 dx O 2ω 1 ρ (w 1 • q 1 ) -∂ ω D(ω 1 )∇ 2 w 1 : ∇ 2 q 1 dx . (48)
According to [START_REF] Koutsawa | A generic approach for the solution of nonlinear residual equations. part i: The diamant toolbox[END_REF], the shape derivatives ρ (Ω)(θ) and D (Ω)(θ) read

ρ (Ω)(θ)(x) = -∂ r H (d Ω ) d Ω (θ)(x) ∆ρ, D (Ω)(θ)(x, ω 1 ) = -∂ r H (d Ω ) d Ω (θ)(x) ∆D(ω 1 ),
where d Ω (θ)(x) is given by equation [START_REF] Hamdaoui | Comparison of non-linear eigensolvers for modal analysis of frequency dependent laminated visco-elastic sandwich plates[END_REF].

Obtaining a descent direction from [START_REF] Jouve | Xd3d[END_REF] is not necessarily easy and we would prefer to recover the classical shape derivative structure of a surface integral on the interface. Corollary 4.1. Suppose that the interface is roughly plane (i.e. the principal curvatures can be neglected) and the thickness parameter ε of the diffuse interface is small. Then a good approximation of (46) is

ω 1 (Ω)(θ) ≈ -Γ θ • n ω 2 1 ∆ρ(w 1 • q 1 ) -∆D(ω 1 )∇ 2 w 1 : ∇ 2 q 1 ds O 2ω 1 ρ (w 1 • q 1 ) -∂ ω D(ω 1 )∇ 2 w 1 : ∇ 2 q 1 dx , (49)
where Γ = ∂Ω ∩ int(O).

For the proof consult [START_REF] Allaire | Multi-phase structural optimization via a level set method[END_REF].

Level-set method for topology optimization

From the previous sections we have all the necessary theoretical ingredients to introduce a gradient method for the minimization of an objective function J(Ω) (see in particular 33 for the definition of shape derivative). As stated in [START_REF] Hadamard | Mémoire sur le problème d'analyse relatif à l'équilibre des plaques élastiques encastrées[END_REF], the shape derivative J (Ω)(θ) can be written as a surface integral on Ω where the function V(x) is given in [START_REF] Fedkiw | Level set methods and dynamic implicit surfaces[END_REF] and (49). Supposing V regular enough so it can be naturally extended to R d , a descent direction is found by defining the vector field

θ = -Vn
where n is a natural extension of the normal (which we implicitly assume for the moment that exists). Then we update the shape Ω as

Ω t = (Id + tθ)Ω,
where t > 0 is a small descent step. Formally we obtain

J(Ω t ) = J(Ω) -t ∂Ω V 2 ds + O(t 2 )
which guarantees the decrease of the objective function. We remark that if V turns out to be not regular enough (as it is the case for some objective functions) there are other possible choices for the definition of the descent direction [START_REF] De Gournay | Velocity extension for the level-set method and multiple eigenvalues in shape optimization[END_REF]. First introduced in [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF], the level-set method is a technique for capturing interfaces which are implicitly defined by the zero level-set of an auxiliary function. In particular, this method has been successfully applied to topology optimization problems [START_REF] Sethian | Structural boundary design via level set and immersed interface methods[END_REF][START_REF] Wang | A level set method for structural topology optimization[END_REF][START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF]. Let the bounded domain O ⊂ R d be the working domain in which all admissible shapes Ω are included. In numerical practice, the domain O will be meshed once and for all. We represent the boundary of Ω by means of a level-set function ψ defined over O such that

   ψ(x) = 0 if x ∈ ∂Ω, ψ(x) < 0 if x ∈ Ω, ψ(x) > 0 if x ∈ (O\Ω). (50)
Under the action of a normal vector field -V(t, x)n(x), where n(x) is a natural extension of the normal given by [START_REF] Fedkiw | Level set methods and dynamic implicit surfaces[END_REF] n(x) = ∇ψ |∇ψ| , the shape Ω evolves according to the Hamilton-Jacobi equation

∂ψ ∂t (t, x) -V(t, x)|∇ψ(t, x)| = 0, ∀t, ∀x ∈ O. (51)
Equation ( 51) is posed in the whole reference domain O, and not only on the boundary ∂Ω, if the velocity V is known everywhere (as will be the case in the sequel). Transporting ψ by (51) is analogous to moving the boundary of ∂Ω (the zero level-set of ψ) along the direction Vn.

A common choice of boundary condition for (51) is ∂ψ ∂n = 0 on ∂O.

As mentioned in [START_REF] Allaire | A level-set method for vibration and multiple loads structural optimization[END_REF], this boundary condition is easy to implement since there is no fixed value to assign for ψ at the boundary, and it also allows the solution of (51) to satisfy a maximum principle. More precisely, new inclusions (or holes) in Ω can appear only by advecting the zero level-set of ψ which changes its topology and cannot come from outside the domain O because of spurious negative values created by the boundary conditions.

The numerical solution of ( 51) is computed with a second order explicit upwind scheme [START_REF] Sethian | Level set methods and fast marching methods: Evolving interfaces in computational geometry, fluid mechanics[END_REF] on a Cartesian grid. Since this scheme is explicit in time, the time stepping must satisfy a CFL condition.

Because of the advection process or numerical diffusion, the level-set function may become too flat or too steep leading to large errors either in the location of its zero level set or in the evaluation of its gradient by finite differences. Therefore, it is usual to regularize it periodically by solving the following problem

∂ψ ∂t (t, x) + sign(ψ 0 )(|∇ψ(t, x)| -1) = 0, ∀t, ∀x ∈ O ψ(t = 0, x) = ψ 0 (x) ∀x ∈ O, (52)
which admits as a stationary solution the signed distance to the initial interface {ψ 0 (x) = 0} [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF].

6. Numerical analysis

Material properties

Within this section we explain how to In order to extend the state equation ( 5) to the whole domain O, we use the same "ersatz material" approach as [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF]. This approach amounts to filling the holes O\Ω by a weak phase mimicking void but avoiding the singularity of the rigidity matrix. This approach in needed in our 3D framework but not in the 2D since this is rather a multi-phase problem. We define an elasticity tensor A (x), which is a mixture of A in Ω and of the weak material mimicking holes in O\Ω, as

A (x) = χ A (x)A, χ A (x) = 1, if x ∈ Ω, δ A , if x ∈ O\Ω. ( 53
)
We also need to apply the same procedure for the material density by introducing a mixture density

ρ (x) = χ ρ (x)ρ, χ ρ (x) = 1 if x ∈ Ω, δ ρ if x ∈ O\Ω. (54)
For eigenfrequency optimization making a correct choice for the threshold parameters δ A and δ ρ is always delicate since a bad combination can yield spurious eigenmodes localized in the ersatz material [START_REF] Allaire | A level-set method for vibration and multiple loads structural optimization[END_REF].

In the case of the composite sandwich structure equation ( 7), the tensor D and the density ρ are evaluated according to the smooth multi-phase approximation [START_REF] Koutsawa | A generic approach for the solution of nonlinear residual equations. part i: The diamant toolbox[END_REF] considering d Ω = ψ, i.e. the current level-set function. Indeed the level-set function ψ is periodically reinitialized through (52). We set ε = 1.5∆x in [START_REF] Schreiber | Nonlinear eigenvalue problems: Newton-type methods and nonlinear Rayleigh functionals[END_REF] where ∆x is the characteristic size of the computational mesh. We remark that in the multi-phase framework we do not need an ersatz material since the whole domain O is covered by the base plate material which makes the rigidity matrix of the composite sandwich always non-singular, no matter the shape of the damping layer.

Solving the discrete non-linear eigenvalue problem

Denote as

T h (ω)u h = 0, u h ∈ V h (55)
the matrix representation of the discretization of problem ( 16) that stems from the weak (or variational) formulation of (5) (the space V h represents the finite element approximation space of (36)). Numerous numerical methods such as the modal strain energy method [START_REF] Johnson | Finite element prediction of damping in structures with constrained viscoelastic layers[END_REF] and the asymptotic numerical method [START_REF] Daya | A numerical method for nonlinear eigenvalue problems application to vibrations of viscoelastic structures[END_REF] have been developed to solve the resulting nonlinear problem. An overview of the different methods of resolution can be found in [START_REF] Koutsawa | A generic approach for the solution of nonlinear residual equations. part i: The diamant toolbox[END_REF][START_REF] Schreiber | Nonlinear eigenvalue problems: Newton-type methods and nonlinear Rayleigh functionals[END_REF][START_REF] Hamdaoui | Comparison of non-linear eigensolvers for modal analysis of frequency dependent laminated visco-elastic sandwich plates[END_REF]. In our case we simply opt for applying the Newton's method to the extended system

F z (ω, u h ) = T h (ω) • u h z • u h -1 = 0, (56)
where z ∈ V h is an arbitrary vector, fixed once and for all, such that z = 1 and z • u h (ω * ) = 0 with u h (ω * ) the exact eigenvector of (55) for eigenvalue ω * . Hence the second equation in (56) represents a normalization condition on u h . The numerical solution (ω 1 , u 1 ) obtained by FEM discretization of the problem ( 15) is chosen as the initial value (ω 0 , u 0 h ) and z = u 0 h . The Newton equation of (56)

F z (ω n , u n h ) + ∂F z (ω n , u n h ) ω n+1 -ω n u n+1 h -u n h = 0
gives the following update rules:

ω n+1 = ω n - 1 z • T -1 h (ω n ) • T h (ω n ) • u n h , u n+1 h = (ω n -ω n+1 )T -1 h (ω n ) • T h (ω n ) • u n h
, where the matrices T -1 h (ω n ) and T h (ω n ) respectively stand for the inverse and the derivative with respect to ω of the FEM matrix T h (ω n ). We iterate until |ω n+1ω n |/|ω n | < tol with tol << 1.

Test cases

Numerical test cases involving the structure eigenvalues in topology optimization are usually not wellposed since less structure implies a higher eigenfrequency as a result of the ersatz material. Hence the optimizer converges to a trivial solution without any material. Typical remedies to solve this classical issue are including non-structural masses, defining the problem as a reinforcement problem or imposing a mass equality constraint (consult for instance [START_REF] Bendsøe | Topology optimization: theory, methods and applications[END_REF]). For the following test cases we apply the first two solutions.

The 3D and 2D eigenvalue computations are performed using Freefem++ [START_REF] Hecht | New development in Freefem++[END_REF] and the 3D results are rendered with XD3D [START_REF] Jouve | Xd3d[END_REF].

The chosen viscoelastic properties for the 3D and 2D cases satisfy the holomorphic condition stated in Section 3 (except on a finite number of poles). That being said, the optimization process remains valid for any other viscoelastic material for which we can assume that the spectrum is discrete. Otherwise, as explained earlier, numerical difficulties, such as spurious eigenvalues, arise when dealing with a continuous spectrum or an accumulation point in the spectrum [START_REF] Dauge | Numerical approximation of the spectra of non-compact operators arising in buckling problems[END_REF]. Remark 6.1. Given a certain operator S : V → V (where V is an adequate function space) and a sequence of discretized operators S n : V n → V n with S n → S in the operator norm when n → ∞ (V n is for instance a finite element subspace of V ), we say that the eigenvalues ω n are spurious in the sense that

ω n ∈ Spectrum(S n ), ω n → ω ∞ /
∈ Spectrum(S).

3D viscoelastic structure

We optimize a three-dimensional cantilever. The working domain O is of size 1m × 2m × 1m (discretized with a 20×40×20 mesh). The viscoelastic material corresponds to 3M ISD112 with complex Young modulus

E c (ω) = (1 + ν c )G 0 1 + 3 k=1 ∆ k ω ω -i z k , G 0 = 0.5MPa, (59)
density ρ V c = 1600 kg/m 3 and Poisson coefficient ν c = 0.5 (consult [START_REF] Xu | Mixed integer multi-objective optimization of composite structures with frequencydependent interleaved viscoelastic damping layers[END_REF] for the values of the not listed constants z k and ∆ k ). A zero displacement boundary condition is imposed on the left side of the domain and four cubic cells on the middle of the right side (heavy tip mass with a material density 100 times heavier than the viscoelastic material) are not subject to optimization (see Figures 4 and5). As mentioned before, this is one of the possible solutions to avoid a trivial optimal shape made only of ersatz material. Indeed, since the ersatz material has by definition a very low density (thus mass), the weight of the structure is minimized while keeping the same eigenfrequency by puling out all the material from the design space. The heavy tip mass is fixed during the optimization process so if the algorithm tried to pull out all the rest of the material from the box, the structure would break (the heavy mass must be connected to the fixed boundary condition). The ersatz material is characterized by the same Poisson's ratio as the viscoelastic material, a smaller Young's modulus by a factor δ A = 10 -2 and a smaller density by a factor δ ρ = 10 -4 (see ( 53) and (54)). The objective function is given as a linear combination of the negative loss factor (since we minimize) and the volume of the structure

J(Ω) = -η + |Ω| |O| , (60)
where = 0.01 is a fixed Lagrange multiplier for the volume constraint. The heavy tip mass is positioned on the rectangular (and not the square) face of the design domain in order to avoid the symmetric bending modes on the X and Y axis. This condition is not necessarily sufficient to enforce the first eigenvalue ω 1 to be simple (our computations are based on this assumption) so we check at each iteration the existence of multiple eigenvectors. The modal loss factors η of the initial and optimal shapes are respectively 0.025 and 0.041. Remark 6.2. In this case Ω corresponds to a uniform isotropic viscoelastic structure such that

A(ω) = f (ω)A 0 ,
where f is a scalar function and A 0 is a real isotropic material independent of ω (both defined accordingly to equation ( 59)). Hence solving the problem ( 5) is equivalent to solving the equation:

div(A 0 e(u)) = λρ u in Ω, (61) which has a countably infinite family of solutions (λ k , u k ) k≥1 , labeled by increasing order of the eigenfrequency. The eigenmodes u k coincide with the eigenmodes of the original viscoelastic problem and we can easily check that p k = u k , so the optimization problem is self-adjoint. In order to compute the modal loss factor η, we determine the original eigenvalue ω 1 as the solution of the non-linear scalar equation: density ρ V c = 1524 kg/m 3 and Poisson coefficient ν c = 0.5. The four parameters a 0 = 332.2, a 1 = 2485.2, c 1 = 0.12, β = 0.47 and the shift factor α(T ) are available in [START_REF] Jones | Handbook of viscoelastic vibration damping[END_REF]. In our case we will consider the temperature T to be constant and equal to 27 • C.

The objective function J(Ω) is the same as (60) but this time with = 10 -1 . The modal loss factors of the initial and optimized shapes are respectively 1.5 × 10 -2 and 2.31 × 10 -2 (see Figure 7). We remark that the optimized shape of the viscoelastic treatment coincides with the one obtained in [START_REF] Kim | Topology design optimization for vibration reduction: reducible design variable method[END_REF] using the SIMP method. The result may not be completely intuitive at first since a significant portion of the damping reinforcement material lays outside the eigenmode largest deflection zone (Figure 8). However, as explained in Section 2.2, the sandwich energy is mainly dissipated as a result of in-plane compression and extension of the damping material under the flexural stress of the base plate. Hence the optimal damping reinforcement material layout and the area with the largest extensional strain energy density overlap (see Figure 9).

Figure 1 :Figure 2 :

 12 Figure 1: Unconstrained-layer damping treatment.

Figure 3 :

 3 Figure 3: Slice of the 3d plate. The shape O belongs to the xy-plane, being x the outward axis in the figure.

Corollary 3 . 1 .

 31 Suppose that the complex elastic tensors A(ω) and D(ω) in (5) and (7) respectively, are holomorphic with respect to ω for ω ∈ D ⊂ C bounded. Suppose also that 0 ∈ D and A(0), D(0) are hermitian positive definite, i.e. A(0) = A(0) H (being the symbol H the conjugate transpose) and ξ H A(0)ξ > 0, ∀ξ ∈ C d×d , where d = 3 (equivalently D(0) with d = 2). Then problems (5) and (7) possess only a finite number of eigenvalues in each compact subset of D.

Lemma 4 . 1 .

 41 Let Ω be a smooth bounded open set and φ

Definition 4 . 4 .

 44 Let Ω ⊂ R d be a Lipschitz bounded open set.

ω 2 -

 2 f (ω)λ 1 = 0, via Newton's method with the initial value ω 0 1 = λ

Figure 4 :

 4 Figure 4: Initial (with a cavity inside) and optimized shapes of a 3D cantilever.

Figure 5 :

 5 Figure 5: The optimized shape possesses a small inner cavity. If the hole is filled with the viscoelastic material, the value of η remains almost constant (actually it slightly diminishes) so the value of the composite objective function (60) increases.The authors verified that the inner cavity is not present when the whole working domain O is used as initial shape in the optimization process.

Figure 6 :

 6 Figure 6: Convergence history for the optimization of the 3D cantilever.

Figure 7 :

 7 Figure 7: Initial and optimized shapes of the composite plate. The aluminum phase is shown in gray and the (superposed) viscoelastic one in black.

Figure 8 :

 8 Figure 8: Real part of the eigenmode w 1 for the final shape.

22

 22

Figure 9 :

 9 Figure 9: In-plane extensional strain energy density (J/m 2 ).

Figure 10 :

 10 Figure 10: Convergence history for the optimization of the 2D composite sandwich plate.

Acknowledgments

This research work was carried out in the framework of IRT SystemX, Paris-Saclay, France, and therefore supported with public funds within the scope of the French Program "Investissements d'Avenir". The authors also thank the financial support of the industrial partners of the TOP project: Safran, Renault, Airbus Group and ESI Group. We would like additionally to show our gratitude to Julien Cortial who read carefully this article providing a very useful feedback and Prof. El Mostafa Daya for fruitful discussions and kind assistance on viscoelastic sandwich modeling.

Computing a descent direction

The gradient of the loss factor

, can be derived from the formulas given in Sections 4.2 and 4.3. Applying the product rule

.

(57)

3D viscoelastic structure

For a given ω 1 , define the constant

) : e(p 1 ). Write [START_REF] Fedkiw | Level set methods and dynamic implicit surfaces[END_REF] as

Then since θ and the normal n are real-valued, (57) reads

The above expression provides directly a normal descent direction θ = V n, with

)γ(ω 1 ) .

Composite sandwich structure

Following the same development, with the definitions

)(ω 1 )∇ 2 w 1 : ∇ 2 q 1 , the same shape derivative (58) applies.


Baptiste Magnier
email: baptiste.magnier@mines-ales.fr

Hasan Abdulrahman
email: hasan.abdulrahman@mines-ales.fr

A Study of Measures for Contour-based Recognition and Localization of Known Objects in Digital Images

Keywords: Distance transform, similarity measures, contours

Usually, the most important structures in an image are extracted by an edge detector. Once extracted edges are binarized, they represent the shape boundary information of an object. For the edge-based localization/matching process, the differences between a reference edge map and a candidate image are quantified by computing a performance measure. This study investigates supervised contour measures for determining the degree to which an object shape differs from a desired position. Therefore, several distance measures are evaluated for different shape alterations: translation, rotation and scale change. Experiments on both synthetic and real images exhibit which measures are accurate enough for an object pose or matching estimation, useful for robot task as to refine the object pose.

I. INTRODUCTION AND CONTEXT

Object detection and recognition are used in many computer vision applications. In the past, several methods have been developed to achieve this in digital images, involving, for example, directly shape context [START_REF] Belongie | Shape context: A new descriptor for shape matching and object recognition[END_REF], patches [START_REF] Crivellaro | Robust 3d object tracking from monocular images using stable parts[END_REF], points of interest [START_REF] Schmid | Evaluation of interest point detectors[END_REF] or region based-methods [START_REF] Zhang | Review of shape representation and description techniques[END_REF]. The shape representation of an object is particularly useful for accurate industrial inspection tasks, where this shape must be aligned with a reference model of the interested object [START_REF] Steger | Machine vision algorithms and applications[END_REF]. A contour-based representation remains a class of methods and exploits only the shape boundary information. In this paper, we focus on distance measures between the acquired features (contours [START_REF] Canny | A computational approach to edge detection[END_REF]) in a candidate image and an ideal contour map model. This represents a supervised evaluation of the shape representation. As the shape of the sought object is already known (by a reference pose, for example), a learning stage is not necessary to estimate the object pose. Consequently, supervised edge detection evaluations compute a score between a ground truth edge map (G t ) and a candidate image (D c ) to achieve this task. Originally, they are used to assess edge detection methods [START_REF] Chabrier | Comparative study of contour detection evaluation criteria based on dissimilarity measures[END_REF], [START_REF] Lopez-Molina | Quantitative error measures for edge detection[END_REF], [START_REF] Magnier | Edge detection: A review of dissimilarity evaluations and a proposed normalized measure[END_REF], [START_REF] Magnier | A Review of Supervised Edge Detection Evaluation Methods and an Objective Comparison of Filtering Gradient Computations Using Hysteresis Thresholds[END_REF]; here, a contour-based localization evaluation of these measures is investigated. Thus, they compute a score of coherence which qualifies the correct object pose possibility.

II. SUPERVISED EDGE DETECTION EVALUATIONS

Several alterations may interfere and disturb the object pose estimation, as translation, rotation or scale change of the interest object. Both their own shape(s) and their contours are altered, as illustrated in Fig. 1. Supervised contour measures are presented below. Then, an evaluation process is fulfilled, determining the degree to which an object shape differs from a desired position as a function of different alterations.

Various evaluation methods have been proposed in the literature to assess different shapes of edges using pixel-based ground truth (see reviews in [START_REF] Dubuisson | A modified Hausdorff distance for object matching[END_REF], [START_REF] Chabrier | Comparative study of contour detection evaluation criteria based on dissimilarity measures[END_REF], [START_REF] Lopez-Molina | Quantitative error measures for edge detection[END_REF], [START_REF] Magnier | Edge detection: A review of dissimilarity evaluations and a proposed normalized measure[END_REF], [START_REF] Magnier | A Review of Supervised Edge Detection Evaluation Methods and an Objective Comparison of Filtering Gradient Computations Using Hysteresis Thresholds[END_REF]). Indeed, a supervised evaluation criterion computes a dissimilarity measure between a ground truth (G t ) and a detected contour map (D c ) of an original image I. In this paper, the closer to 0 the score of the evaluation is, the more the segmentation is qualified as suitable. To assess an edge detector, the confusion matrix remains a cornerstone in boundary detection evaluation methods. Comparing pixel per pixel G t and D c , the 1st criterion to be assessed is the common presence of edge/nonedge points. A basic evaluation is composed of statistics by combining G This type of assessment is useful for region segmentation evaluation [START_REF] Crum | Generalized overlap measures for evaluation and validation in medical image analysis[END_REF], but, a reference-based edge map quality measure requires that a displaced edge should be penalized in function of FPs and/or FNs and of the distance from the position where it should be located [START_REF] Magnier | Edge detection: A review of dissimilarity evaluations and a proposed normalized measure[END_REF] [START_REF] Magnier | A Review of Supervised Edge Detection Evaluation Methods and an Objective Comparison of Filtering Gradient Computations Using Hysteresis Thresholds[END_REF], as shown in Fig. 2.

F oM (G t , D c ) = 1 - 1 max (|G t | , |D c |) • p∈Dc 1 1 + κ • d 2 Gt (p) κ ∈ ]0; 1] F oM revisited [15] F (G t , D c ) = 1 - 1 |G t | + β • F P • p∈Gt 1 1 + κ • d 2 Dc (p) κ ∈ ]0; 1] and β ∈ R + Combination of F oM and statistics [17] d 4 (G t , D c ) = 1 2 • (T P -max (|G t | , |D c |)) 2 + F N 2 + F P 2 (max (|G t | , |D c |)) 2 + F oM 2 (G t , D c ) κ ∈ ]0; 1] Edge map quality mea- sure [18] D p (G t , D c ) = 1/2 |I|-|Gt| • p∈F P 1 - 1 1 + κ•d 2 Gt (p) + 1/2 |Gt| • p∈F N 1 - 1 1 + κ•d 2 T P (p) κ ∈ ]0; 1] Symmetric FoM [9] SF oM (G t , D c ) = 1 2 • F oM (G t , D c ) + 1 2 • F oM (D c , G t ) κ ∈ ]0; 1] Maximum FoM [9] MF oM (G t , D c ) = max (F oM (G t , D c ) , F oM (D c , G t )) κ ∈ ]0; 1] Yasnoff measure [19] Υ (G t , D c ) = 100 |I| • p∈Dc d 2 Gt (p) None Hausdorff distance [20] H (G t , D c ) = max max p∈Dc (d Gt (p)), max p∈Gt (d Dc (p)) None Maximum distance [12] f 2 d 6 (G t , D c ) = max   1 |D c | • p∈Dc d Gt (p), 1 |G t | • p∈Gt d Dc (p)   None Distance to G t [21][12][8] D k (G t , D c ) = 1 |D c | • k p∈Dc d k Gt (p), k = 1 for [21] and [12] k ∈ R + Oversegmentation mea- sure [22] Θ (G t , D c ) = 1 F P • p∈Dc dG t (p) δTH k
for [START_REF] Odet | Scalable discrepancy measures for segmentation evaluation[END_REF]:

k ∈ R + and δ T H ∈ R * + Undersegmentation measure [22] Ω (G t , D c ) = 1 F N • p∈Gt dD c (p) δTH k
for [START_REF] Odet | Scalable discrepancy measures for segmentation evaluation[END_REF]:

k ∈ R + and δ T H ∈ R * + Relative Distance Error [12], [23], [24] RDE k (G t , D c ) = k 1 |D c | • p∈Dc d k Gt (p) + k 1 |G t | • p∈Gt d k Dc (p), k ∈ R + , k = 1 for [12], k = 2 for [23], [24] Symmetric distance [12][8] S k (G t , D c ) = k p∈Dc d k Gt (p)) + p∈Gt d k Dc (p) |D c ∪ G t | , k = 1 for [12] k ∈ R + Baddeley's Delta Metric [25] ∆ k (G t , D c ) = k 1 |I| • p∈I |w(d Gt (p)) -w(d Dc (p))| k k ∈ R + and a convex function w : R → R Magnier et al. measure [26] Γ(G t , D c ) = F P +F N |Gt| 2 • p∈Dc d 2 Gt (p) None Complete distance mea- sure [9] Ψ(G t , D c ) = F P +F N |Gt| 2 • p∈Gt d 2 Dc (p) + p∈Dc d 2 Gt (p) None λ measure [28] λ(G t , D c ) = F P +F N |Gt| 2 • p∈Dc d 2 Gt (p) + min |G t | 2 , |Gt| 2 T P 2 • p∈Gt d 2 Dc (p) None Edge Mismatch Measure (EM M ) [27] EM M (G t , D c ) = 1 - T P T P + ω • p∈F N δ Dc (p) + • p∈F P δ Gt (p) M dist ∈ R + , D max ∈ R + , ω ∈ R + , ∈ R + . δ Dc (p) = d Dc (p), if d Dc (p) < M dist D max , otherwise and δ Gt (p) = d Gt (p), if d Gt (p) < M dist D max , otherwise.
In [START_REF] Sezgin | Survey over image thresholding techniques and quantitative performance evaluation[END_REF]:

M dist = 0.025•|I|, D max = |I|/10, ω = 10/|I|, = 2.
Table 1 reviews the most relevant measures involving distances. Thus, for a pixel p belonging to the candidate contour D c , d Gt (p) represents the minimal Euclidian distance between p and G t . These types of distance measures play an important role in image matching and may be used to determine the degree of resemblance between two objects [START_REF] Huttenlocher | A multi-resolution technique for comparing images using the hausdorff distance[END_REF] [START_REF] Dubuisson | A modified Hausdorff distance for object matching[END_REF]. To achieve this, if p belongs to G t , d Dc (p) corresponds to the minimal distance between p and D c , Fig. 2 illustrates the difference between d Gt (p) and d Dc (p). Mathematically, denoting (x p , y p ) and (x t , y t ) the pixel coordinates of two points p and t respectively, thus d Gt (p) and d Dc (p) are described by:

     for p ∈Dc: d G t (p)= Inf (xp -xt) 2 + (yp -yt) 2 , t ∈ Gt , for p ∈Gt: d Dc (p)= Inf (xp -xt) 2 + (yp -yt) 2 , t ∈ Dc .
These distance functions refer to the Euclidean distance.

On the one hand, some distance measures are specified in the evaluation of over-segmentation (i.e., distances of FPs), for example: Υ, D k , Θ and Γ; see also [START_REF] Dubuisson | A modified Hausdorff distance for object matching[END_REF], [START_REF] Magnier | Edge detection: A review of dissimilarity evaluations and a proposed normalized measure[END_REF], [START_REF] Magnier | A Review of Supervised Edge Detection Evaluation Methods and an Objective Comparison of Filtering Gradient Computations Using Hysteresis Thresholds[END_REF]. On the other hand, the Ω measure assesses an edge detection by computing only under-segmentation (distances of FNs). Other edge detection evaluation measures consider both distances of FPs and FNs [START_REF] Abdulrahman | From contours to ground truth: How to evaluate edge detectors by filtering[END_REF]. A perfect segmentation using an oversegmentation measure could be an image including no edge points. As demonstrated in [START_REF] Magnier | A Review of Supervised Edge Detection Evaluation Methods and an Objective Comparison of Filtering Gradient Computations Using Hysteresis Thresholds[END_REF], another limitation of only over-and under-segmentation evaluations are that several binary images can produce the same result. Therefore, a complete and optimum edge detection evaluation measure should combine assessments of both over-and under-segmentation, as H, ∆ k , f 2 d 6 , S k , RDE k , Ψ and λ.

One of the most popular descriptors is named the Figure of Merit (F oM). This distance measure has an advantage because it ranges from 0 to 1, where 0 corresponds to a perfect segmentation [START_REF] Abdou | Quantitative design and evaluation of enhancement/thresholding edge detectors[END_REF]. Nonetheless, for F oM, the distance of the FNs is not recorded and contours having a small displacement compared to their desired positions are strongly penalized as statistic measures (detailed in [START_REF] Magnier | A Review of Supervised Edge Detection Evaluation Methods and an Objective Comparison of Filtering Gradient Computations Using Hysteresis Thresholds[END_REF]). Several evaluation measures are derived from F oM: F , d 4 , EM M , MF oM, SF oM and D p . Contrary to F oM, the F measure computes the distances of FNs but not of the FPs. In addition, the d 4 measure depends particularly on T P , F P , F N and ≈1/4 on F oM, but d 4 penalizes FNs like the F oM measure. Note that EM M computes a score different from 1 if there exists at least one TP (cf. experiments). Otherwise, SF oM and MF oM take into account both distances of FNs and FPs, so they can compute a global evaluation of a contour image. However, MF oM does not consider FPs and FNs at the same time, contrary to SF oM . Another way to compute a global measure is represented by D p . Nevertheless, D p is more sensitive to FNs than FPs because of its huge coefficient (cf. Table 1).

A second measure widely computed in matching techniques is represented by the Hausdorff distance H, measuring the mismatch of two sets of points [START_REF] Huttenlocher | A multi-resolution technique for comparing images using the hausdorff distance[END_REF]. This measure is useful in object recognition, the algorithm aims to minimize H, reporting the mismatch of two shapes [START_REF] Hemery | Comparative study of localization metrics for the evaluation of image interpretation systems[END_REF] [START_REF] Paumard | Robust comparison of binary images[END_REF]. It is well known that this distance could be strongly deviated by only one pixel positioned sufficiently far from the pattern, there are several enhancements of H presented in [START_REF] Huttenlocher | A multi-resolution technique for comparing images using the hausdorff distance[END_REF][12] [START_REF] Baudrier | Binary-image comparison with local-dissimilarity quantification[END_REF]. As pointed out in [START_REF] Dubuisson | A modified Hausdorff distance for object matching[END_REF], an average distance from the edge pixels in D c to those in G t is more appropriate, like S k , RDE k or Ψ. Moreover, S k takes small values in the presence of low level of outliers, whereas the score becomes large as the level of mistaken points increases [START_REF] Dubuisson | A modified Hausdorff distance for object matching[END_REF][8] but is sensitive to remote misplaced points [START_REF] Magnier | Edge detection: A review of dissimilarity evaluations and a proposed normalized measure[END_REF]. Otherwise, derived from H, the Delta Metric (∆ k ) [START_REF] Baddeley | An error metric for binary images[END_REF] intends to estimate the dissimilarity between each element of two binary images, but is highly sensitive to distances of misplaced points. The λ measure penalizes highly FNs compared to FPs, because the more FNs are present in D c , the more D c the desirable object becomes unrecognizable.

Finally, an objective assessment of these measures (and others) has been performed by varying the thresholds on contours of images obtained by filtering techniques, see [START_REF] Chabrier | Comparative study of contour detection evaluation criteria based on dissimilarity measures[END_REF][9][10], [START_REF] Abdulrahman | A New Objective Supervised Edge Detection Assessment using Hysteresis Thresholds[END_REF] [START_REF] Magnier | A Review of Supervised Edge Detection Evaluation Methods and an Objective Comparison of Filtering Gradient Computations Using Hysteresis Thresholds[END_REF]. Theoretically, the minimum score of the measure corresponds to the best edge map, compared to the ground truth. Here, a study is led as a function of the position of the candidate object contours which may be corrupted by undesirable pixels caused by object translation, rotation of scale change (or/and due to noise and/or blur in real images).

III. EVALUATIONS AND EXPERIMENTAL RESULTS

Object localization, recognition and matching are practical tools in the computer vision domain. Usually, they are designed to meet the requirements of industrial alignment applications (automated optical inspection, industrial automation, target search, visual servoing...). The distance measures detailed above are compared here with the specific purpose of contour-based object localization evaluation. Such measures are tested on synthetic and real image in order to assess if they could be useful and accurate enough to recognize and match objects. The two subsections above present both experiments on synthetic and real images. 22 error measure are tested:

Dice * , F oM, F , d 4 , MF oM, SF oM , D p , Υ, H, f 2 d 6 , D k k=2 , Θ, Ω, RDE k=1 , RDE k=2 , S k k=1 , S k k=2 , ∆ k , Γ, Ψ,
λ and EM M , showcasing the advantages and drawbacks of each one. The tests carried out in the experiments are intended to be as complete as possible, and thus as close as possible to reality. Note that the Matlab code of the distance measures as F oM, D k , S k and ∆ k are available at http://kermitimagetoolkit.net/library/code/. Others are available on the mathworks website: https: //fr.mathworks.com/matlabcentral/fileexchange/ 63326-objective-supervised-edge-detectionevaluation-by-varying-thresholds-of-the-thinedges. Firstly, to verify if the measure has the required properties, different alterations are applied to create synthetic localization results simulating real ones.

A. Localization and Recognition of Synthetic Objects

In order to quantify the reliability of a measure of dissimilarity, an edge map of a synthetic shape is affected by various alterations: rotation, translation and scale change. Thus, it allows to verify if the evolution of the provided score obtained by a measure is in line with the expected behavior: small error for close shapes and large penalty for shapes that differ. Note that there exist matching techniques imposing grid of specific window on the image, typically in order to search small objects in the image [START_REF] Huttenlocher | A multi-resolution technique for comparing images using the hausdorff distance[END_REF][32] [START_REF] Boughnim | Hand posture classification by means of a new contour signature[END_REF] or line tackers [START_REF] Deriche | Tracking line segments[END_REF].

Here the distance measures consider both full images of the same size: the ground truth (with an isolated object) and the candidate object contour. For the quantification of similarity, special consideration should be given to the ability of a skill measure to penalize important dissimilarities and ignore the unimportant ones; as an example few spurious isolated undesirable pixels must not disturb a measure if the object is close to its desired location. To do that, a synthetic shape is created (see Fig. 3(a)). Candidate shapes are also created to simulate movements, thereupon compared with the desired location (ground truth). For assessment, a validation measure is sought so that it can appropriately penalize differences in the three continuous deformations: 1) rotation, 2) translation and 3) scale. Thus, a curve is obtained for each measure and each alteration type. Only some curves are reported in this paper, curves which are not reported are similar to other ones.

1) Rotation:

The test is performed by rotating the control shape incrementally up to a total rotation of 2π, as illustrated in Fig. 3. Then, the metric distances are computed between the rotated contour and the reference shape. The shape of the curve of the measure scores is expected to be roughly symmetric around π. It is noteworthy that locations of local extrema of the curve depend on the shape of the object. Usually, a minimum is computed for distance measures at 180 • , due the rough symmetry of the shape, as for H. For example, in Fig. 3(d), some contours of the rotated shape are aligned with the desired location, creating two minima at 160 • and 200 • for the λ measure. So, this measure penalizes strongly FNs, so it computes an error at 180 • larger than before or after, because various details of the shape are missing. However, the Dice measure computes a pixel per pixel score and does not record distances, so its behavior is unstable. Also, F oM measure and its derivatives are also unstable and very sensitive to small rotations, minima are not clearly detectable. 2) Translation: In this test, the contour shape is gradually translated by moving the shape away from its original location along a horizontal straight line, cf. Fig. 4 (the results are the same for translations in other directions). The resulting measure scores have a minimum at 0, where the two contours are collocated. Scores obtained for measures as H and RDE grows linearly, as a function of the shape displacement. Concerning Dice or EM M , the obtained curves have obvious discontinuities, demonstrating the limited sensitivity of these measures to translation. EM M obtains a score of 1 when no TP exists. Also, F oM and its derivatives are sensitive to displacements, with a jump after few pixel displacements and with score attaining its maximum (value 1) when no common points between the two shapes. The over segmentation of Υ is not monotonous when not TP are present. Other measures behave desirably, with a good response: increasing monoton-ically. Note that the λ measure decreases at the end of this experiment because there exists less pixels in the candidate shape than in the image containing all the object contours.

3) Scale Changement: For this experiment, the shape undergoes a decrease in size with the maximal scale 8th the original size (see Fig. 5). Scores from the measures are computed between the ground truth and the rescaled shape. For all shapes, excepted Dice, each measure exhibits the expected behavior of being strictly decreasing. However, the graph of EM M has sharp discontinuities that exhibit its unstable response to scaling, because its responses are 1 without TP. The scores of the F oM and its derivatives are not decreasing enough, excepted where the two contours are collocated. Finally, λ behaves correctly converging throw 0 until a scale of 3, and continue to converge less rapidly after.

B. Localization and Recognition of Real Objects

Experiments on color real images were also performed. In Figs. 6 and8, thin edges are extracted using the Canny edge detector (σ = 1) [START_REF] Canny | A computational approach to edge detection[END_REF][36] and a non-maximum suppression: the selected pixels are those having gradient magnitude at a local maximum along the gradient direction, which is perpendicular to the edge orientation [START_REF] Ziou | Edge detection techniques: An overview[END_REF]. Here, by moving the camera, the aim is to determine when the object is at the desired position in the image using thin binary edges as features, the scores must converge to 0 as summarized in Fig. 7. This desired position corresponds to the object in the last video frame, as shown in Figs. 6(c) and 8(c), each frame may be corrupted by numerous FPs and the candidate object may contain FNs. The first video contains 264 frames, whereas the second video is composed of 116 frames of size 720×1280. The ground truth corresponds to binary boundaries of the desired position of the object, represented in blue pixels in Figs. 6(a) and 8(a). Green pixels represent TPs, red points are FPs whereas blue pixels which are also G t are tied to FNs. These features are dilated with a structural element of size 3×3 for a better visualization.

For the first video, object contours are easily extracted, with spurious undesirable points at the end of the video, as illustrated in Fig. 6(a)-(c). Consequently, only edge points detected out of the candidate object shape may disturb the contour-based localization. Indeed, curves tied to H, ∆ k , Υ and S k behave stochastically (see the caption of Fig. 6 concerning other measures). Also, Dice and F oM scores converge only to lower values when the candidate object is close to the desired location. Only Ω and λ measures have the intended behavior for this sequence, as summarized in Tab.2.

The second video is heavily corrupted by a random noise on the each color plane (SNR≈11dB). These disturbances create spurious pixels in the edge detection process, but especially, the edges of the candidate object are not well localized or not present. Therefore, the majority of measures behaves stochastically, and the last frame does not represent the best score for them (see Fig. 8 and its caption). For the F oM measure, usually its scores decrease, even though the curve converges considerably only for the last frames. On the EM M measure converges rapidly, but, remains constant after few number of frames. The Ω measure is relatively constant, excepted for the last frames where it decreases. Lastly, the λ measure behaves as expected with a minimum at the end.

IV. CONCLUSION

This study presents an assessment of measures for contourbased recognition and localization of known objects. Thus, the Dice and 20 distance measures are evaluated trough different shape alterations: translation, rotation and scale change. Experiments on real images exhibit which measures are accurate enough for an object pose or shape matching estimation. Clearly, a measure involving false negative distances remains more accurate than other techniques, as the under-segmentation measure Ω. However, when some parts of the candidate image are missing, but detected close to their desired positions, they are not taken into account by Ω. Alternatively, the λ measure remains well defined for this assessment. Missing edges are penalized higher than spurious Table 2. Reliability of the reviewed edge detection evaluation measures.

Meas.

Tr. ones because some isolated points disturb the shape location, as the majority of the measures, cf. Tab.2. Finally, λ does not need any tuning parameter or specific window and may be useful to validate a visual servoing process.

Rot. Sc. Real Dice * ✗ ✗ ✗ ✗ F oM ✗ ✗ ✗ ✗ F ✗ ✗ ✗ ✗ d4 ✗ ✗ ✗ ✗ SF oM ✗ ✗ ✗ ✗ MF oM ✗ ✗ ✗ ✗ Dp ✗ ✗ ✗ ✗ Υ ≈ ✓ ✓ ✗ H ✓ ✓ ✓ ✗ f2d6 ✓ ✓ ✓ ✗ Meas. Tr. Rot. Sc. Real D k ≈ ✓ ✓ ✗ Θ ✓ ✓ ✓ ✗ Ω ✓ ✓ ✓ ≈ RDE ✓ ✓ ✓ ✗ S k ✓ ✓ ✓ ✗ ∆ k ✓ ✓ ✓ ✗ Γ ≈ ✓ ✓ ✗ Ψ ≈ ✓ ✓ ✗ λ ✓ ✓ ✓ ✓ EMM ≈ ✓ ✗ ✗
The studied measures are tied to edge binarized images (edge detection methods are not discussed, see [START_REF] Magnier | Edge detection: A review of dissimilarity evaluations and a proposed normalized measure[END_REF], [START_REF] Abdulrahman | From contours to ground truth: How to evaluate edge detectors by filtering[END_REF], [START_REF] Magnier | A Review of Supervised Edge Detection Evaluation Methods and an Objective Comparison of Filtering Gradient Computations Using Hysteresis Thresholds[END_REF]). In practice, illumination changes modify seriously the image contrast and, additionally to the noise and the blur, may be an obstacle for shape recognition. To overcome this problem, the combination of edge-based measures and nonlinear illumination changes measure [START_REF] Steger | Similarity measures for occlusion, clutter, and illumination invariant object recognition[END_REF] Frame number

Fig. 1 .Fig. 2 .

 12 Fig. 1. Example of translation, rotation and scale change of Dc.

  t and D c . Afterwards, denoting | • | as the cardinality of a set (e.g. |G t | represents the number of edge pixels in G t ), all points are divided into four sets (cf. Figs.1): • True Positive points (TPs): T P = |G t ∩ D c |, • False Positive points (FPs): F P = |¬G t ∩ D c |, • False Negative points (FNs): F N = |G t ∩ ¬D c |, • True Negative points (TNs): T N = |¬G t ∩ ¬D c |. Various edge detection evaluations involving confusion matrices have been developed , cf. [8][9][11]. The Dice measure [13] is one example: Dice * (G t , D c ) = 1 -2•T P 2•T P +F N+F P .

( a )Fig. 3 .Fig. 4 .Fig. 5 .

 a345 Fig. 3. Examples of localization metrics behaviors for a rotation alteration for each 5 • : from 0 • to 360 • . The legend is available in Fig. 1(d). F oM, F , SF oM and MF oM obtain almost the same scores. Score evolutions for are the same for Dice * , d 4 and Dp. Also Υ behaves like D k , Γ, Ψ and Θ which is the perfect symmetry of Ω. Scores tied to S k k={1,2} , RDE k={1,2} , f 2 d 6 , ∆ k and Ψ are also similar. EM M behaves lire λ.

Fig. 6 .

 6 Fig. 6. Localization metrics behaviors on real experiment. The image in (b) represents the edge movements as a function of the time (from blue to red). F oM, F , SF oM and MF oM obtain almost the same scores which are close to Dice * , d 4 and Dp evolution. Scores tied to H, ∆ k , Θ, Γ, Ψ, Υ, D k , S k k=2 , RDE k=1 , RDE k=2 and f 2 d 6 are similar. Score evolutions are the same for S k k=1 , EM M and λ.

Fig. 7 .

 7 Fig. 7. Expected behavior of a measure scores.

Fig. 8 .

 8 Fig. 8. Localization metrics behaviors on real experiment with noisy images. The image in (b) represents the edge movements as a function of the time (from blue to red), showing the huge number of noise pixels along the video. F oM, F , SF oM and MF oM obtain almost the same scores which are close to Dice * , d 4 and Dp evolution. The measure ∆ k behaves like the Hausdorff distance, stochastically along the video without convergence. Also Υ behaves like D k , Θ, S k k=2 , RDE k=1 , RDE k=2 and f 2 d 6 . Scores tied to Ψ, Γ , EM M and S k k=1 are also similar.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Localization metrics behaviors on real experiment tied to Fig. 6 in the original paper.

Table 1 .

 1 List of error measures involving distances, generally: k = 1 or k = 2 , and, κ = 0.1 or κ = 1/9 .

	Error measure name	Formulation	Parameters
	Pratt's FoM [14]

  could be deeply investigate. This document presents additional results on real images.The showed curves are tied to Fig.6and Fig.8in the original paper.

		1					1					1							1
	A Study of Measures for Contour-based Recognition and Localization of Known Objects in Digital Images: Additional Results 0.4 0.5 0.6 0.7 0.8 0.9 Measure score 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Measure score 0.4 0.5 0.6 0.7 0.8 0.9 Measure score 0.4 0.5 0.6 0.7 0.8 0.9 Measure score
	0.3					0.2					0.3							0.3
	0.2					0.1					0.2							0.2
	Measure score	0.9 1 0.3 0.4 0.5 0.6 0.7 0.8	Baptiste Magnier and Hasan Abdulrahman IMT Mines Alès, LGI2P, 6. avenue de Clavières 30100 Alès, France 100 150 200 250 Frame number 50 100 150 200 250 Frame number 50 100 150 200 250 Frame number (a) Dice * (b) F oM (c) d4 0.95 1 0.45 0.5 Abstract. 50 0.65 0.7 0.75 0.8 0.85 0.9 Measure score 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Measure score	Measure score	0.25 0.05 0.1 0.15 0.2	50	100 Frame number 150 200 ( d ) F	250
		0.2					0.6					0.05
		0.1					0.55					0							0
			50	100	150	200	250	50	100	150	200	250	50	100	150	200	250			50	100	150	200	250
				Frame number			Frame number			Frame number					Frame number
			(e) SF oM			(f) MF oM				(g) Dp					(h) EM M
		550
	Measure score	250 300 350 400 450 500
		200
		150
			50	100	150	200	250	50	100	150	200	250
				Frame number