Catherine Quilliet

Numerical deflation of beach balls with various Poisson's ratios: from sphere to bowl's shape

We present a numerical study of the shape taken by a spherical elastic surface when the volume it encloses is decreased. For the range of 2D parameters where such a surface may model a thin shell of an isotropic elastic material, the mode of deformation that develops a single depression is investigated in detail. It occurs via buckling from sphere toward an axisymmetric dimple, followed by a second buckling where the depression loses its axisymmetry through folding along portions of meridians. For the thinnest shells, a direct transition from the spherical conformation to the folded one can be observed. We could exhibit unifying master curves for the relative volume variation at which first and second buckling occur, and clarify the role of the Poisson's ratio. In the folded conformation, the number of folds and inner pressure are investigated, allowing us to infer shell features from mere observation and/or knowledge of external constraints.

Introduction

Let's consider a thin shell of an elastic isotropic material, such as a beach ball, and deflate it. What would be its shape ? This question is not restricted to garrulous familial shores : fundamental and applied physics nowadays presents legions of easily deformable soft objects, and knowing what governs their shapes gives the powerful possibility of inferring mechanical properties from simple observations, without contact. Among these deformable objects, an increasing number derives from spherical symmetry, that is omnipresent at scales where surface effects overcome volume forces such as gravity. The numerical study presented in this paper discusses the shapes taken by spherical thin shells of isotropic materials when their inner volume is decreased by a significant amount. Such a systematic and quantitative study will help deciphering conformations observed in e.g. Soft Matter (lock-and-key colloids [1], multiwall capsules [2], particles design through evaporation [3,4]), galenics (encapsulation [5]), microfluidics (microtanks [6]) or medicine (ultrasound contrast agents [7]), under the action of an external pressure or other possibly isotropic fields such as concentration in evaporation/dissolution phenomena, which also shapes objects in Nature [8].

When an elastic spherical shell has its inner volume lowered, it first deforms through in-plane compression that respects the spherical symmetry. Then it undergoes a symmetry breaking in order to relax a high stretch energy into much lower bending energy, by reversion of a spherical cap (creation of an axisymmetric depression, or "dimple"). The onset of this sudden transition, or buckling, under external pressure was studied long ago by Pogorelov and Landau [START_REF] Pogorelov | Bending of surfaces and stability of shells[END_REF][START_REF] Landau | Theory of elasticity[END_REF]. They showed that the dimple should nucleate over a critical outside/inside pressure difference scaling like Y 3D d R 2 (in what follows, we will refer to this latter quantity by ∆P Landau ), where d is the shell thickness, R its radius and Y 3D the Young modulus of the material that makes it up. Its edge (or "rim") has a transversal extension √ dR. One of the key of their 1 calculation being the assumption that buckling occurs for dimples such that maximum deflection is of order d, spherical geometry imposes then that √ dR is also the radius of the dimples that forms. Besides, classical buckling analysis provided dependence of the buckling pressure with the Poisson's ratio [9,[START_REF] Pogorelov | Bending of surfaces and stability of shells[END_REF][START_REF] Landau | Theory of elasticity[END_REF].

Results focusing on deformations through further deflation are, mainly, more recent. A few months ago, stability analysis allowed a detailed study of buckling toward axisymmetric conformations [START_REF] Knoche | [END_REF]. Experimental [13,3,14,15,16,1,17] and numerical [3,18,17] deflation studies showed shapes holding several dimples, also called "multiple indentation". These conformations compete with experimental observations of shapes holding a single depression [14,16], possibly losing axisymmetry [20] or exhibiting folding perpendicularly to the rim [21,15]. Similar shapes are observed in shells under a point load [22,23] or pressed against a wall [24,23]. Secondary buckling by folding of the single depression, also called "polygonal indentation", under isotropic constraint was numerically retrieved with surface models [21,18]. Thin shells with a single depression, either axisymmetric or polygonal, appear to present a conformation of lower energy than the metastable multiple indentation [21,18,17]. In the case of an axisymmetric dimple, this can be easily understood since elastic energy mainly concentrates in dimple edges as bending energy, with an energy per edge length that weakly varies with dimple size. Hence dimples coalescence lowers the total elastic energy [START_REF] Quilliet | Calculations of[END_REF]. Nevertheless, more than one dimple may nucleate if the deflation is rapid enough, leading to metastable [3,18,17] multi-indented shapes. The term "rapid" is to be taken on a wide acception here. Experimentally, it may correspond to situations where dissipation (due to material viscosity or to fluid flows accompanying the deformation) prevents dimple growth, which favors secondary nucleation once ∆P Landau is reached, and where subsequent kinetics prevents thermally activated coalescence between adjacent dimples. Numerically, minimization may reproduce such metastable situations [3,18,17], since (i) large volume increments favors the creation of extra dimples, by making difficult to find the cooperative displacement of vertices that corresponds to rim rolling in dimple growth (ii) depending on the way curvatures are calculated, energy barriers that prevents from dimple coalescence may be overcome or not. For "slow" deflations, a single dimple can appear and grow, or freshly nucleated dimples may coalesce into a single one. Such "slow" deflation provides an axisymmetric bowllike shape, that may undergo under further deflation a transition toward a non-axisymmetric depression, i.e. polygonal indentation [21,18].

We present here a systematic numerical study of such "slow" deflations leading to shapes with a single depression. In this purpose, we used a surface model taking into account recent developments, presented in section 2. We clearly expose the correspondence between 2D parameters of the model surface, and 3D properties of the real object of nonzero thickness, expliciting the role of the different significative parameters. Particular emphasis is put on a parameter often underconsidered: the Poisson's ratio.

The whole study allows to determine parameters of importance for the transitions sphere → axisymetric bowl (section 3) and axisymmetric bowl → polygonal indentation (section 4), both for the detailed shape in polygonal indentation, and for inner pressure. Furthermore, we took particular care to provide empirical dependence laws for practical use.

Surface model

Surface model, where out-of-plane and in-plane deformations are formally uncoupled, is for long considered as valid to describe the deformation of thin sheets (plates or shells) [START_REF] Landau | Theory of elasticity[END_REF][START_REF] Audoly | Elasticity and Geometry: From hair curls to the nonlinear response of shells[END_REF]. For thin sheets without spontaneous curvature (i.e. an elementary surface portion of the sheet, freed from constraints exerted by surrouding material, remains flat at equilibrium), the energies per surface unit that are to be considered in this surface model are of two kinds: firstly, a curvature term that can express 1 2 κc 2 + κg [START_REF] Helfrich | [END_REF][START_REF] Landau | Theory of elasticity[END_REF], where c = 1 R1 + 1 R2 and g = 1 R1 × 1 R2 are respectively the mean and Gaussian curvatures (R 1 and R 2 being the local principal curvature radii), and κ and κ are respectively the mean and Gaussian curvature constants [START_REF]There may be discrepancies between different communities about the definitions of both κ and c. Here we used definitions that rule in the Soft Matter community, and more particularly in physics of lipids vesicles: c = 1 R1 + 1 R2 , where R 1 and R 2 are the main algebraic curvature radii[END_REF]. The other term may be written, in a Hookean linear model: 1 2 ǫ ij K ijkl ǫ kl , where ǫ ij and K ijkl respectively represent the two-dimensional strain and elasticity tensors for in-plane deformations.

For an homogeneous and isotropic surface, the nonzero terms of the two-dimensional elasticity tensor are

K xxxx = K yyyy = Y2D 1-ν 2 2D , K xxyy = K yyxx = ν2D Y2D 1-ν 2 2D
and K xyxy = K yxyx = Y2D 1+ν2D , with Y 2D the two-dimensional Young modulus and ν 2D the two-dimensional Poisson ratio, which is comprized between -1 and 1 [START_REF] Landau | Theory of elasticity[END_REF]. This in-plane elasticity term can be rewritten as Y2D 2(1+ν2D ) Tr(ǫ 2 ) + ν2D (Trǫ) 2

1-ν2D

for the sake of concision.

In a linear approximation, the relation between the 2D parameters and the 3D features of the plate (Young modulus Y 3D , Poisson's ratio ν 3D , thickness d) with zero boundary tangential constraints [START_REF] Landau | Theory of elasticity[END_REF] is expressed as (detailed e.g. in [24] or [7]):

ν 2D = ν 3D = ν (1)
Y 2D = Y 3D d (2)
κ = Y 3D 12 (1 -ν 2 ) d 3 (3)
κ = (ν -1) κ = - Y 3D 12 (1 + ν) d 3 (4)
Since for bulk materials the maximum value of ν 3D is 1 2 for thermodynamic reasons [START_REF] Landau | Theory of elasticity[END_REF], one can notice that the range of 2D Poisson's ratio that effectively describes a thin plate of an isotropic material is limited to a maximum value of 1 2 . In other terms, even a thin plate of an incompressible isotropic material cannot behave as an incompressible surface (where ν 2D = 1), thanks to the possibility of having its thickness varied. On the other limit, Poisson's ratio can reach -1 as a lower value, but negative values correspond to less common "auxetic" materials.

Conversely the thickness of the plate, as a function of 2D parameters, writes:

d = 12 (1 -ν 2 ) κ Y 2D (5)
For describing surfaces with asymmetric properties, the notion of "spontaneous curvature" was introduced by W. Helfrich [START_REF] Helfrich | [END_REF]. It was recently shown that to describe the deformations of an initially stress-free thin shell of radius R, the three contributions (in-plane, mean curvature and gaussian energy, that may have a non-vanishing part even for closed surfaces, depending on the definition of a nonzero spontaneous curvature) can be rewritten in an easily computable way as [START_REF] Marmottant | [END_REF]:

E elastic = cst + shell surf ace 1 2 κ (c -c * 0 ) 2 + 1 2 ǫ ij K ijkl ǫ kl + γ ef f dS (6)
with c * 0 = 1+ν R being the effective spontaneous curvature, and

γ ef f = (1-ν 2 )κ 2R 2
an effective surface tension.

This expression is slightly different from the one used in [21], hence we will quantitatively discuss, in the results, modifications induced by the use of this more complete expression.

More generally, we will consider the influence of sphere size through the use of the adimensionalized Föppl-von Kármán number [START_REF] Lidmar | [END_REF]

γ = Y2D R 2 κ
, that gives the order of magnitude of the ratio between inplane and out-of plane deformation energies [24]. An elastic surface with the energy given in equation ( 6) can effectively describe a thin shell of an isotropic material if 12 1 -ν 2 /γ ≪ 1, in addition to the condition ν 2D ≤ 1 2 . In this range, γ roughly scales like R d 2 . Out of this range, such a surface model does not correspond to any thin shell of an isotropic material ; it can nevertheless describe different types of objects, e.g. gel-phase vesicles [17] that can hence be considered as thin shells of non isotropic materials.

Numerical experiments are performed by minimizing the elastic energy as expressed in equation ( 6) for different inner volumes, with the free software Surface Evolver [31]. A whole in-silico deflation experiment (from V init = 4 3 πR 3 to ≈ 0.2 × V init ) is realized through a succession of different equilibrium states, these latter found according the process described in [21], and calculated successively for inner volumes decreased by steps of at maximum 2% of the initial volume (steps amplitude is reduced in some situations in order to avoid the nucleation of secondary dimples).

3 First-order transition toward axisymmetric depression Deflation of a spherical elastic surface, at imposed either volume or external pressure, causes an abrupt buckling from the spherical conformation in order to release in-plane compressive stress. The purpose here is to compare the numerical approach described in the previous paragraph to known features of this buckling, for the range of parameters that scans the generality of thin shells of isotropic material.

Buckling pressure

First buckling relaxes in-plane constraints, and causes drastic drop of the inside/outside pressure difference ∆P = P ext -P int . Figure 1 displays typical evolutions for the pressure difference: first a linear increase followed by a drop at first buckling, after which pressure difference varies in a much lesser extent. Linear behaviour is expected before the first buckling due to the relation between pressure and elastic energy P ext -P int = ∂E elastic /∂ (∆V ) (detailed in [7]), and quadratic dependence of E elastic with ∆V = V init -V (see e.g. [21]).

In the simulations presented here, designed not to be stuck in multi-indentation conformations of higher energy, or other less stable, the first buckling leads to a single axisymmetric dimple. It is expected to happen when the external overpressure ∆P = P ext -P int reaches the critical value [9,[START_REF] Knoche | [END_REF]:

∆P c = 2 3 1 -ν 2 -1/2 × Y 3D d R 2 (7)
As expected, first buckling in our simulations effectively occurs at a pressure difference of order

∆P Landau = Y 3D d R
2 (section 1). For a given Poisson's ratio, the incertitude due to discrete volume increments does not allow to conclude that ∆P buckling1 /∆P Landau is affected by γ (see e.g. figure 1 displaying several cases at ν = -0.5 ). The effect of the Poisson's ratio is displayed on figure 2: it shows that the pressure which induces buckling in our simulations quantitatively follows the theoretical equation 7, which reinforces the validity of our approach.

Buckling volume

The mechanism of buckling in axisymmetric conformations was quite recently investigated in detail by Knoche et al [START_REF] Knoche | [END_REF], with the study of various metastability branches. At first significative order, their calculations show that for the trivial isotropic ("spherical") deformation, the relative volume variation due to an external overpressure ∆P = P ext -P int expresses:

∆V V = 3 (1 -ν) 2 × R d × ∆P Y 3D
Hence the deflation at buckling pressure writes [9,[START_REF] Knoche | [END_REF]:

∆V V buck 1 = 3 1 -ν 1 + ν × d R (8)
Figure 3 presents values of ∆V V buck 1 from numerical simulations, as a function of a combination of γ and ν translated in 3D parameters: here also the theoretical equation 8 is quantitatively retrieved. One may notice (since by essence ∆V V buck 1 ≤ 1) that the sphere→bowl transition vanishes for d R ≥ (1 + ν) /3 (1 -ν), which indicates a destabilization of the axisymmetric bowl for the most auxetic materials. This was qualitatively expected from the κ divergence in this limit, which makes curvature deformations prohibitive compared to in-plane compressions [START_REF]A possible expression of the 2D compressibility modulus is χ 2D = Y2D 2(1-ν) = Y3D ×d 2(1-ν)[END_REF].

Using the 2D parameters of the surface model, relation 8 also writes:

∆V V buck 1 = 6 (1 -ν) γ -1 2 (9)
Extending our purpose out of the range of 2D parameters that effectively can describe a thin shell of isotropic material, we may remark that extrapolation to ν = 1 induces vanishing of ∆V V buck 1 , and hence destabilization of the spherical conformation : this limit corresponds to incompressible surfaces, that will necessarily undergo a deformation implying curvature even for the smallest volume decreases, since area variation is prohibited. For the thickest shells, i.e. γ 80, shape transition occurs not any more through sudden inversion of a spherical cap, but by slowly deforming into an ovoid, that flattens at the location of future depression under further deflation. We did not, here, specifically study this extreme behaviour.

For practical purposes, we may notice that the large range of parameters explored shows that in terms of volume, the onset of buckling mainly depends on the relative thickness d R , with only a weak influence of the Poisson's ratio. Since this latter ranges between 0 and 1 2 for common materials (i.e. non auxetic, with ν ≥ 0), the prefactor of d R varies between √ 3 and 1, which is much narrower than the range in d R that can be explored.

4 Second-order transition toward polygonal depression

Location of the transition

In the axisymmetric bowl shape, global bending of the rim on the equator costs in-plane deformation: extension on the outer side of the rim, and compression on the inner one. For the thinnest shells, compressive stress parallel to the equator leads to a secondary buckling, where the inner side of the rim undulates to adapt to axial compression (fig. 4, left), forming folds, or "wrinkles", that deform the axisymmetric depression into a roughly polygonal shape (fig. 4, right). Such a conformation mainly involves curvature deformations, much less energetic [24] than compression energy that quadratically increases with ∆V V . Fig. 5 shows how elastic energy dispatches between in-plane and out-of-plane deformation energies in a typical numerical deflation. Wrinkles match with the rim through a zone of high curvature that has a folding role similar to what realizes the apex of d-cones [START_REF] Pauchard | [END_REF]34,35], except that the surface is not developable but is spherical, hence the concept of "s-cones" proposed by Reis and Lazarus [36].

Secondary buckling from axisymmetric bowl shape to polygonal indentation, quite smooth, is harder to detect than the first one (fig. 5). The corresponding relative volume variation ∆V V buck 2 is determined on one hand by the maximum deflation before loss of axisymmetry, and on the other hand by the ∆V V at which the rim presents convex zones under axial observation (as shown on lower part of Fig. 5, subfigure c). Figure 6 presents a typical shape phase diagram for ν = -0.5, with three distinct zones: spherical coformation, axisymmetric bowl and wrinkled depression. One may notice that for the thinnest shells, the incremented deflation we numerically performed shows a direct transition from the sphere to the wrinkled bowl. Relative volume variation at first buckling nevertheless obeys a single power-law on the whole range of relative thicknesses. We focus in this section on the second buckling, hence considering only deflations where wrinkles appear on an already existing axisymmetric conformation. 6), i.e. 1 2 ǫ ij K ijkl ǫ kl . Grey continuous line: elastic energy linked to the first significative term :

1 2 κ (c -c * 0 ) 2 .
This term does not vanish at ∆V V = 0 because c * 0 = c 0 ; in the total elastic energy, it is counterbalanced by the cst and the effective surface term expressed by γ ef f in equation ( 6). Effective surface energy (not represented here) varies at maximum by -0.33 κ , and on an amplitude 0.03 κ in the nonspherical conformations. Grey zone indicates the second transition from axisymmetric bowl shape to depression with inner wrinkles ("polygonal indentation" [18]), determined as explained in the text. First transition from sphere to axisymmetric bowl occurs through abrupt decrease of the in-plane deformation energy, at ∆V V = 0.03. logarithms [START_REF]Plotting ln ∆V V buck 2 versus ln d R , as for the first buckling, excessively particularizes points at ν = -0.8[END_REF] provides with a correlation coefficient of -0.99:

∆V V buck 2 = 8470 × γ -1.085 (10)
This expression is compatible with previous results obtained for ν = 1/3 with slightly different numerical models: the shell without gaussian curvature evoked in section 2 [21], and the spring model without spontaneous mean curvature of reference [18] (both are represented on figure 7 for their range of validity).

Extrapolating equation 10 up to ∆V V = 1 suggests that this secondary buckling does not happen, i.e. single indentation keeps its axisymmetry, below a threshold value γ c,buck2 = 4170. In tridimensional parameters, relation [START_REF] Pogorelov | Bending of surfaces and stability of shells[END_REF] expresses:

∆V V buck 2 = 571 × d/R √ 1 -ν 2 2.17 (11)
The dependence in ν for non-auxetic materials is even weaker than for the axisymmetric buckling since it plays at maximum by a factor 4/3. Besides, for the range of Poisson's ratio studied here, there is no axisymmetric conformation to be expected for d/R < 0.003. Relation (11) also implies that wrinkles are not expected when d/R 0.054 × √ 1 -ν 2 : a particular consequence is that in wrinkles prevention, a very auxetic material (with ν → -1) may help.

Characterization of buckled shapes with the number of wrinkles

For the thinnest shells that undergo polygonal indentation, the most conspicuous feature is the number W of wrinkles, or s-cones. Figure 8 shows the evolution of W for a typical deflation: first does W decrease while the freshly nucleated and still very flat depression hollows and enlarges, then it increases again. Data are quite scattered: there is a typical noise of order ±1 on W , that has no observable correspondence in smooth energy curves (fig. 7). In order to decrease data scattering, we calculated W def lated as the average value of W between ∆V V = 0.53 and ∆V V = 0.76 (these values have been choosen in order to cover, for all the simulations, a significant range of relative volume variation before autocontact, this latter happening around ∆V V ≈ 0.9). Values of W def lated are comparable with results from the previous model, which did not take gaussian curvature into account, and indicates a scaling law in γ -1/4 [21,7]. The dependence of W def lated with γ and ν is shown on figure 9, expressed in 3D parameters. It shows a scaling in d R -1/2 , which provides clues on the typical transversal size l of the s-cones (presented on figure 4, left). Since s-cones stand alongside one another on a length which is of the order of an equator, we can estimate l as 2πR W def lated . Hence best fit of W def lated with dependence in d R -1/2 (figure 9) can be expressed as:

l ≈ 6.7 √ dR ( 12
)
This result is fully comparable to the wrinkles wavelength 4.7 √ dR that can be calculated from recent results by Vella et al [19] on the indentation of strongly pressurized shells. As shown in reference [START_REF] Landau | Theory of elasticity[END_REF], √ dR arises naturally from balancing the bending and in-plane deformation energy of a small deformation on a spherical shell. Recent results by [START_REF] Knoche | [END_REF] showed that √ dR scales also for the rim width even in large axisymmetric depressions ; we confirm here that it governs also other types of large deformations such as s-cones transversal size in polygonal depressions.

Other postbuckling features

First buckling and its consequences were presented in section 3. For application purpose, we may interest to the inside/outside pressure difference after its drop at first buckling. Careful examination of the numerical data revealed supplementary features, uncorrelated with the occurrence of the secondary buckling previously exposed : two types of behaviour clearly appear in the evolution of the reduced pressure ∆P ∆P Landau during deflation. For d R 0.014, pressure difference ∆P presents the type of evolution calculated by [START_REF] Knoche | [END_REF], i.e. quasiplateauing after buckling (variation of about 15% during the whole deflation, plus some occasional dispersion due to numerical procedure), up to autocontact. Furthermore, an order relation is respected: at every volume step, the ratio ∆P ∆P Landau weakly increases when γ decreases. This is to be observed in fig. 1, for the 4 curves corresponding to the highest relative thicknesses.

For the thinnest shells ( d R 0.012 in the simulations performed), ∆P ∆P Landau regularly re-increases with deflation after the pressure drop, crossing successively the curves at smaller γ's (as shown on the 2 "thinnest" curves of fig. 1).

In order to extract a general behaviour from these differents observations, we focused on ∆Pmin ∆P Landau , the minimum value of ∆P ∆P Landau after buckling. Figure 10 shows that the cross-over between the two regimes around d R ≈ 0.014 also corresponds, for each ν, to the minimum of relative pressure drop after the first buckling. In the plateauing regime, ∆Pmin ∆P Landau stands for the plateauing value for d R 0.014 ; it shows a power-law of the type ∆Pmin ∆P Landau = a (ν) × d R 0.5 . Similarly to what was done in section 3, and since a (ν) appears to be even (curves at ν = 0.5 and ν = -0.5 almost mix up on figure 10), we looked for a prefactor of the form 1 -ν 2 x , minimizing x for the best fits at d R > 0.014. This led us to propose the master curve presented on figure 11, of formula ∆Pmin Y3D × 1 -ν 2 0.773 = 0.75 × d R 2.5 . We do not have for the moment theoretical clues to justify these two successive fitting operations, but (i) it allows to describe numerical results in a very condensed way for d R ≥ 0.014 , and (ii) for all the shells these reduced values impressively gather on a single curve, for γ ranging from 8 × 10 2 to 4.7 × 10 5 , and for ν between

Conclusion

Systematic numerical study of the buckling of a spherical shell, in the conformation with a single depression, allows to sketch the influence of the different geometrical or elastic parameters through quite simple theoretical or phenomenological laws. The surface model can be translated in 3D parameters, that are the shell's thickness, and the two elastic parameters of the material that compose it: Young modulus and Poisson's ratio.

At imposed volume, the Young modulus does not play on the shape. Results showed that the first transition (toward axysimmetrically buckled shape), and the second one, with appearance of wrinkles, or "s-cones", is mainly driven by d R for non-auxetic (i.e. with positive Poisson's ratio) materials. For auxetic materials, the Poisson's ratio may have a determining importance, by strongly displacing transitions toward higher values of the relative volume variation, up to possible vanishing. Decreasing the Poisson's ratio down to very negative values stabilizes spherical deflation at the expense of dimples creation, and axisymmetric dimples against appearance of wrinkles.

The number of wrinkles indicate a dependence in d R -1/2 , that confirms √ dR as the accurate scaling for elastic deformations of elastic spherical surfaces. The Young modulus scales pressure features: critical inside/outside pressure difference that triggers first buckling, and plateauing pressure after buckling. Detailed behaviour, that is shown to reduce to a master curve, opens the possibility for two different wrinkling regimes.

Figure 1 :

 1 Figure 1: Pressure difference ∆P = P ext -P int adimensionalized by ∆P Landau = Y 3D d R 2 (or 12 Y 2D κ 1 -ν 2 1/2 /R 2 in 2D parameters), as a function of the relative volume variation, for surfaces of similar Poisson's ratio ν = -0.5. Dotted green: γ = 1.17 × 10 5 ( d R = 8.8 × 10 -3 ); blue: γ = 6.07 × 10 4 ( d R = 1.22 × 10 -2 ), interrupted orange: γ = 3.22 × 10 4 ( d R = 1.67 × 10 -2 ), magenta: γ = 1.17 × 10 4( d R = 2.31×10 -2 ), light blue: γ = 9.33×10 3 ( d R = 3.11×10 -2 ), parma: γ = 4.67×10 3 ( d R = 4.39×10 -2). Notice scale switch at ∆V V = 0.2. Points (corresponding each to a minimization) emphasizes two curves with typical behaviour after buckling: increasing (green curve), and plateauing (orange) .

Figure 2 :

 2 Figure 2: Pressure (∆P ) buck 1 = (P ext -P int ) buck 1 at first sphere→bowl buckling, adimensionized by ∆P Landau . Error bars represent the range of this quantity for γ between 4.67 × 10 3 and 2.33 × 10 5 . Line: theoretical value 2 3 1 -ν 2 -1/2 , from equation 7.

Figure 3 :

 3 Figure 3: Relative volume variation ∆V V buck 1 at the sphere→bowl buckling (lower and upper limit of the error bars are indicated respectively by upward and downmard triangles). White: ν = 0.5 ; light gray: ν = 0 ; dark grey: ν = -0.5 ; black: ν = -0.8. Grey line: ∆V V buck 1 = 6 (1 -ν) γ -1 2 (this is equation 8 expressed in 2D parameters). Dotted line: relation established in [21] for ν = 1/3: taking into account the gaussian curvature elastic energy increases volume threshold by at maximum 20% for the thinnest shells of the shells previously studied. Illustrative inserts display deflated spherical surfaces at γ = 4666 and ν = -0.8 (hence d R = 0.0304), with relative volume variations respectively ∆V V = 0.161 (spherical) and 0.167 (buckled, section view, same scale).

Figure 4 :

 4 Figure 4: Left: schematisation of a spherical surface near complete deflation. Inner part of the rim endures a compressive stress in the direction indicated by the red double arrow. Relaxation of in-plane deformation occurs, for the thinnest shells, via undulation deformation, generating wrinkles (folds) all along the rim (examples schematized with blue interrupted lines). Length ℓ stands for lateral extension of the wrinkles. Right: simulation with γ = 9.33 × 10 3 , ν = 0.5 ( d R = 0.031) and ∆V V = 0.562 ). Maxima of undulation are stressed with blue interrupted lines ; circle locates the zone of high curvature that forms the apex of the s-cone.

  Figure 7 displays how ∆V V buck 2 varies with the Föppl-von Karmán number γ for different Poisson's ratios. Data indicate a dependence in γ, but scattering prevents from concluding on an influence by ν ; linear regression on

Figure 5 :

 5 Figure 5: Deformation energies at different relative volume variations (volume step ≈ 0.6% of the initial volume), adimensionalized by the curvature constant κ, for γ = 9.33 × 10 3 and ν = 0.5 ( d R = 0.031). Interrupted line: Hookean in-plane deformation energy, linked to the second significative term of eq. (6), i.e.1 2 ǫ ij K ijkl ǫ kl . Grey continuous line: elastic energy linked to the first significative term :1 2 κ (cc * 0 ) 2 .

Figure 6 :

 6 Figure 6: Shapes in the relative thickness / relative volume variation phase diagram, for: ν = -0.5. Illustrations: shells with d/R = 0.031 and ν = -0.5 ( ∆V V =0, 0.20 and 0.55).

Figure 7 :

 7 Figure 7: Points: relative volume variation at which the axisymmetric depression becomes "polygonal", versus 12/γ (which is to be identified with d/R when ν = 0). White: ν = 0.5 ; light grey: ν = 0 ; dark grey: ν = -0.5 ; black: ν = -0.8. Wide grey line: equation 11. Continuous lines indicate the location of the first-order buckling for different Poisson's ratio (from equation8), under which only the spherical conformation is to be found. Non-continuous lines: in their domain of validity, equations describing the axisymmetric/wrinkled bowl transition published in previous works for ν = 1/3 (see text): dotted black refers to[21], interrupted grey to[18].

Figure 8 :

 8 Figure 8: Number of wrinkles (s-cones) held by the single depression after the secondary buckling of figure 5. Poisson's ratio ν = 0.5; Föppl-von Karman number γ = 6.06 × 10 4 ( d R = 0.0122). Lower part: conformations at points indicated on the main figure.

Figure 9 :

 9 Figure 9: Number W def lated of wrinkles (s-cones) at the end of a numerical deflation (averaged between ∆V V = 0.53 and ∆V V = 0.76). error bars are taken as the standard deviation on this range, with a minimum value of ±0.5. White: ν = 0.5 ; light grey: ν = 0 ; dark grey: ν = -0.5 ; black: ν = -0.8. Continuous line: W def lated = 0.940 × d R

Figure 10 :

 10 Figure 10: Minimum inside/outside pressure difference after buckling (cf figure 1), adimensionalized by Landau pressure, logarithmic representation. Black: ν = 0.5; dark grey: ν = 0; light grey: ν = -0.5; white: ν = -0.8. Interrupted line separates the two types of evolution of the pressure after the first buckling: increasing or plateauing (see figure1 ). Black line indicate slope 0.5.

Figure 11 :

 11 Figure 11: Pressure master curve after first buckling: ∆Pmin Y3D × 1 -ν 2 0.773 versus d R . Black squares: ν = 0.5; black diamonds: ν = 0; light grey diamonds: ν = -0.5; white diamonds: ν = -0.8. Continuous line: ∆Pmin Y3D × 1 -ν 2 0.773 = 0.75 × d R 2.5 .

Acknowledgments

The author thanks K. Brakke for developing and maintaining the Surface Evolver software, including invaluable interactions during this work, and P. Marmottant and F. Quéméneur for fruitful discussions.


Hichem Arioui
email: arioui@iup.univ-evry.fr

SAID MAMMAR Abderrahmane Kheddar

Said Mammar

STABLE SHARED VIRTUAL ENVIRONMENT HAPTIC INTERACTION UNDER TIME-VARYING DELAY

Keywords: Virtual environment haptics, Varying time delayed control, stability and robustness

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Virtual reality techniques, refer typically to humanin-the-loop or human centered advanced simulation or prototyping systems. The original feature of the concept lies in the multi-modality of the man-machine interaction, which involves all human sensory modalities. Among these capabilities, the haptic modality is of prime importance when it's a matter to allow the human operator to experience honest manipulation and touching of virtual objects with realistic sensations of stiffness, roughness, temperature, shape, weight, contact forces, etc. These physical parameters are collected then interpreted by the human haptic modality through a direct touch and motion of, let say, human hand. Virtual environments are visually rendered to the human operator through screens, head mounted displays and other up-to-date advanced visual interfaces. Headphones are used to display 3D virtual sound. In the contrary to vision and auditory, haptic displays are active. Indeed, to render and display forces, the interfaces must be able to both constraint human desired motions and, to apply forces on the involved human part (e.g. hand). These interfaces are typically robotic devices that are capable: (i) to collect desired human motion or desired human applied force to be sent to the VE engine part of the simulation state update, and (ii) to display, to the human operator, subsequent virtual forces, computed thanks to computer haptics algorithms (collision detection, dynamic contact and reaction force computation, etc.). Applications of force reflection or force feedback are actually spreading to many domains. Among the well known ones: interactive surgical simulators, interactive driving simulators, interactive games, VE based teleoperation. A great demand is also experienced in virtual industrial prototyping. The last issue would extend to concurrent engineering and needs the potentiality to allow haptic interaction among a group of users sharing the same VE over a network. It is well known from physiological and psychological data of the haptic modality and from the haptic control theory that the haptic loop requires a high bandwidth of around 1 kHz to guarantee the stability of the haptic interaction, and more importantly, to make a coherent feedback between the visual and the haptic scenes. Developing a network protocol that can provide sufficient bandwidth with minimum latency to a group of distant users is a challenging problem [START_REF] Buttolo | Architectures for shared haptic virtual environments[END_REF] and physics-based models that simulate haptic interactions among users have begun to be developed [START_REF] Basdogan | Real-time simulation of dynamically deformable finite element models using modal analysis and spectral lanczos decomposition methods[END_REF]. Yet one of the important problems of haptic feedback, even if only one user interacts with the VE engine, is time delay. The difficult nature of some tasks, the lack of knowledge on user abilities and behaviors, the problem of developing a universal controller for stable haptics could be also too complex. Obviously time delay during the transfer and processing of data may easily result in unstable forces and can be harmful to the user.

To the best knowledge of the authors, there is no work addressing the stability of VE delayed force feedback interaction, since in most known applications, the user is not distant from the interactive simulation engine. This papers proposes a simple and efficient solution to deal with this problem.

The linear case with constant time-delay

VE force reflecting techniques borrow a lot to the early teleoperation systems. Up-to-now, time delay is still known to be one of the most sever problem in force reflecting teleoperators. Many solutions have been proposed to deal with this problem. Some of the most attractive ones are based on passivity derived from scattering network theory [START_REF] Anderson | Asymptotic stability for force reflecting teleoperators with time delay[END_REF], [START_REF] Niemeyer | Stable adaptive teleoperation[END_REF], [START_REF] Stramigioli | Geometric Scattering in tele-manipulation of port controlled Hamiltonian systems[END_REF] and [START_REF] Yokokojhi | Bilateral teleoperation uner time-varying communication delay[END_REF]. Others are based on a passive transformation of power parameters (namely velocity and force) into waves. Other techniques use known classical control theory to derive stable controller from Lyapunov criteria [START_REF] Eusebi | Force reflecting telemanipulators with time-delays: Stability analysis and control design[END_REF] and [START_REF] Oboe | A Design and Control Environment for Internet-Based Telerobotics[END_REF], we can state other works as [START_REF] Leung | Bilateral Controller for Teleoperators with Time Delay via µ-Synthesis[END_REF] using µ-synthesis and [START_REF] Otsuka | Bilateral telemanipulation system with time delay based on force-sum driven virtual internal models[END_REF] with notion of virtual time delay, etc. Although Smith prediction method is known since 1957 [START_REF]AC on t ro l le rt oo v e rc o m ede a d time[END_REF], it has not been implemented in the early time-delayed force reflecting teleoperation systems. The reason that prohibits the use of Smith prediction approach lies in the practical impossibility to predict mainly (i) the remote environment behavior and, (ii) the operator desired trajectories, since they are given on-line. Concerning virtual reality applications, since most controllers come from teleoperation experience [START_REF] Miller | Passive implementation for a class of static nonlinear environments in haptic display[END_REF] and [START_REF] Adams | Stable haptic interaction using the excalibur force display[END_REF], it was not surprising to notice that Smith prediction was not investigated as a potential solution for time-delayed VE haptic feedback controllers.

The originality of the proposed solution is in the some-how prediction of the master part within the remote part [START_REF] M A M M A R ,A N Dt .H A M E L | A smith-prediction based haptic feedback controller for time delayed virtual environments systems[END_REF]. Hence, the developed equations lead to a scheme where only the master model appears and also the estimation of the time delay is necessary. The term "somehow prediction" is used to signify that in fact the proposed solution is not really a prediction since only the master model is required, which means that no prediction of operator behavior or trajectory is needed. However, the upwards and forwards time delays must be known.

Without loss of generality and to better understand the concept, a simple LTI model of a VE haptic interface is considered. tual coupling and VE positions, F e is the VE computed force, F h is the operator applied force on the device, C(s) is the commonly used virtual coupling [START_REF] Adams | Stable Haptic Interaction with Virtual Environments[END_REF] and [START_REF] Colgate | Issues in the haptic display of tool use[END_REF], which guarantee unconditional stability of the haptic interaction system in the absence of timedelay, finally, τ 1 and τ 2 are respectively upwards and forwards constant time delays. The closed loop transfer function of the haptic system without the proposed controller is given by:

F e (s) F h (s) = M (s)E(s)e -sτ 1 1+e -s(τ 1 +τ 2 ) E(s)(M(s)+C(s)) (1)
This transfer function has an infinite eigenvalues as the time delay element is present in the characteristic equation. This may consequently imply an instability of the whole haptic interaction.

The proposed solution to overcome this instability is designed within the colored box of figure 1. The controller uses the process model of the haptic display it performs like a local feedback loop within the local remote environment (real or virtual). The resulting transfer function of the global system is a stable haptic feedback with a delayed input F h :

F e (s) F h (s) = M (s)E(s)e -sτ 1 1+E(s)(M (s)+C(s)) (2)
We can notice that, when using the proposed controller there is no more delay items in the characteristic equation of the closed loop system, equation 2. As stated before, the main advantage of this controlprediction scheme is in using the model of the haptic reflecting device only. The latter is well know and its parameters well identified. However, the controller requires to estimate both upwards and forwards timedelay.

Obviously, this example can be easily generalized to any kind of haptic devices and mainly teleoperation systems. In this last case, E(s) represents the linear model of the remote system interacting with its environment.

The linear case with time-varying delay

As previously stated, the proposed method needs:

• a good knowledge of the master model, and

• the estimation of the time delay Estimation of the model does not require complicated techniques. Well known model estimation methods (namely those developed in robotics) can be applied. Estimation of time delay is also easily made, especially when the delay is constant. Simple networking commands (such as a ping command) achieve the matter. Practically, the time delay may fluctuate since neither public or usual networking access protocols nor the computer haptics algorithms 2 (in the frame of VE haptics) can guarantee time-determinism.

A more adequate practical implementation

When we look more the proposed controller (figure 1) we notice that in fact, the desired remote position is corrupted by a local closed loop on the obtained contact force. Indeed, the actual x e is not directly the desired master position x m , but the delayed x m minus the outputs of the local closed controller based on the master model, the delayed master model and the obtained contact force.

C(s) M(s)

e -≥Οs e -≥Νs

F h x m x c - - + + F e E(s) M(s)+C(s)
x e + -

M(s)+C(s)

+ + Figure 2: A practical implementation of the controller.

At the beginning, the controller have been implemented as it is. But from a practical and a simple ob-2 collision detection algorithms and dynamic force computation.

servation, the structure of the controller makes possible an interesting extension which:

• avoids the estimation of time-delay, and

• makes a straightforward extension to timevarying time delay.

• unburden the VE from the buffering in order to compute the controller by sharing the computation on both sites. This is obtained simply as depicted in figure 2. The controllers in the figures 1 and 2 are identical. Nevertheless, the second implementation highlights that it is no more necessary to estimate time delay, and more importantly: the behavior of time delay may have no effect on the stability of the system. The next section gives a more generic discussion and proves the given assertions.

The dynamic model of an haptic display can be approximated in a linear form 3 , considering an apparent mass M and friction B:

M ẍm = F h + τ

A simple controller of the form:

τ(t)=-B ẋ - αF e (t -τ 2 (t)) leads to F h (t) -αF e (t -τ 2 (t)) = M ẍm + B ẋm (3)
where x m , ẋm and ẍm are respectively the Cartesian space position, speed and acceleration, F h and F e denote respectively human and VE forces applied to the haptic device, α is a simple gain taken as 1 in this paper 4 . Time delays uses for the proof of the stability with the proposed control scheme are variable. The fundamental idea is to emulate a passive behavior of the haptic device and the transmission channel. In the second comparator of figure 3 we have x h = x m + x c , using the equality F e (tτ 2 (t)) = M ẍc + B ẋc are obtains:

F e (t -τ 2 (t)) = M [ẍ h -ẍm ]+ B [ ẋh -ẋm ] ( 4
)
Where M is the estimated apparent mass, B is the friction estimate. Equation 4 is used to cancel the effect the delayed control F e (tτ 2 (t)) from the haptic device (master) position x m to be sent to the slave site. Just after the transmission channel τ 1 (t), at the third comparator we have:

F e (t)= M [ẍ h (t -τ 1 (t)) -ẍmr ]+ B [ ẋh (t -τ 1 (t)) -ẋmr ] (5)
where x mr is the after-master position transmitted to the slave site. Equation 4 delayed by τ 1 (t) substracted into equation 5 leads to:

F e (t -τ 1 (t -τ 2 (t))) -F e (t)= M [ẍ mr -ẍm (t -τ 1 (t))] + B [ ẋmr -ẋm (t -τ 1 (t))] (6)
Finally, equation 3 is also delayed by τ 1 (t),a n dt h e obtained force F e (tτ 1 (tτ 2 (t))) substituted in equation 6 leads to:

F h (t -τ 1 (t)) -F e (t)= M ẍmr + B ẋmr +(M -M )ẍ m (t -τ 1 (t)) + (B -B) ẋm (t -τ 1 (t)) (7)
If we assume that the estimation error of the apparent mass and friction is zero:

½ M -M =0 B -B =0
then equation 7 takes the following form:

F h (t -τ 1 (t)) -F e = M ẍmr + B ẋmr (8)
This last equation exhibits a passive behavior of the equivalent new master side. The correction is equivalent to delay the input F h . Assuming that the virtual environment is passive, a fundamental property is that the feedback interconnection of passive systems is again passive [START_REF]L 2 -Gain an Passivity Techniques in Nonlinear Control[END_REF], it ensues from it, that the haptic interaction is stable.

Simulation results

This section presents simulation results of the developed controller. The haptic display is a one DOF actuated arm with apparent mass m =0 .2kg and friction of about b =3 Ns/m. The contact will be performed between the rigid virtual pen and a virtual walls of high stiffness K e . In this first simulation, time delays are taken constant but different, indeed τ 1 =1 s e cand τ 2 =0 .5sec. Figure 4 shows the tracking and force feedback behavior when the operator interacts with a VE of stiffness K e . It is assumed that collision detection and force computation are performed simply and do not cost additional time delay. Figure 4 shows the result obtained of a simulating virtual contact. The operator applies a sinusoidal force F h which drops the master and probe positions to increase until a contact is made between the probe and the wall. This is done when the virtual probe position reaches 20cm. From this time the local VE controller C(s) guarantees the local stability of the virtual interaction and the calculation of the virtual force to feed back to the operator. One can notice that when the contact is made, the mater velocity v m drops to zero and the fed force F e (the controller) increases accordingly to F h during the contact. The position discrepancy, when the contact is made, is unavoidable whatever is the controller or the approach (unless a very prediction is made in the master side), this is due to the undergo physical time-delay. Nevertheless, the virtual probe position x e is stably maintained by the operator during the whole contact time. Several simulations are conducted with multiple hard and viscous contacts, they show that the behavior of the force feedback interaction is stable whatever the time delay. Obviously, one must not suspect that functional performances are acceptable for an actual use in the presence of important time delays.

Figure 5 shows the result of the previously set up system under similar parameters and time-varying delay. The variable time delay τ 1 (t) is represented on the same figure, it satisfies the causality constraints: tτ 1 (t) > 0 and τ 1 (t) < 1 ∀t. In the simulation 0. performed: 1) time-varying τ 1 (t) with τ 2 (t) constant, 2) time-varying τ 2 (t) with τ 1 (t) constant and, 3) timevarying τ 1 (t) and τ (t). Obtained simulation results show that τ 2 (t) dynamics has no effect on the stability of the overall force feedback system. This results make true the derived theoretical assumption stating that the adopted new implementation of the developed controllers cancels the effect of time delay. By analogy to the constant time delay case, since the time delays exponential are not present (after correction) in the characteristic equation. So in case 2, the simulation shows a stable behavior of the system, as if τ 2 (t) was constant i.e. the dynamic of τ 2 (t) is compensated by the controllers. For case 1 and 3, clearly the dynamic of τ 1 (t) affects the dynamic of the force reflecting system. Figure 5 shows clearly that in the transitions between two time delays behavior of τ 1 (t) i.e. varying and constant, the force reflecting response switches between two overall system corresponding behaviors. The transitions seem to be abrupt but do not affect the overall stability of the system. Similarly to the constant time delay linear case, this is related to the fact that the τ 1 exponential is still present in the closed loop transfer function numerator. Indeed the behavior of τ 1 (t) affects the behavior of force feedback but not its stability, as proven in theory, equation 8.

Robutness Analysis

Robustness analysis is performed in the case of constant time delay. The time delays has been approximated using Padé's fourth order transfer function approximation, that is:

exp(-sτ) = lim n→∞ µ 1 -τ 2n s 1+ τ 2n s ¶ n
A root locus of the closed loop transfer function has been performed by varying the controller's parameters that is to say. Figure 6 shows that within a determined margins, force reflection is still stable although some light oscillations appear in the master position x m and the reflected force

Conclusion

The paper presents a master-model based controller designed to stabilize delayed force feedback systems. The proposed method is based on an astute implementation of a somehow Smith prediction scheme, which requires only the haptic device model and do not necessitate the estimation of both (upwards and downwards) delays. Simulation results confirm a stable force reflection from the VE in presence of constant and also time-varying delays. A robustness analysis of the proposed controller has also been conducted. The error margins that guarantee the stability of fed back forces are found to be wide enough to allow using a linear model of the haptic interface based only on a apparent mass and friction estimation.

Comparing to wave-based approaches, the proposed solution is more transparent to the user, since there is no additional corrupting damping as engendered from the transformation of force and flux parameters into waves. The price to be paid is in the position discrepancy between the master and the virtual avatar when the contact is made which may be more important in our controller case comparing to wave-based methods. In fact, in wave-based method, the artificial damping increases with speed (in free motion), which prevents important master-slave position discrepancies, but there is an additional force felt which is not directly related to actual remote contact forces, that is to say more stable but less transparent.

Future work is focused in improving performances in virtual environment haptics. A prediction within the master site is possible based on computer haptics algorithms.

Figure 1 Figure 1 :

 11 Figure 1: The master-model-based Smith prediction principle in the frame of a nominal LTI haptic feedback architecture.

Figure 3 :

 3 Figure 3: The actual implementation scheme.

Figure 4 :

 4 Figure 4: Simulation of delayed virtual contact with a stiff wall and force feedback.

Figure 5 :

 5 Figure4shows the result obtained of a simulating virtual contact. The operator applies a sinusoidal force F h which drops the master and probe positions to increase until a contact is made between the probe and the wall. This is done when the virtual probe position reaches 20cm. From this time the local VE controller C(s) guarantees the local stability of the virtual interaction and the calculation of the virtual force to feed back to the operator. One can notice that when the contact is made, the mater velocity v m drops to zero and the fed force F e (the controller) increases accordingly to F h during the contact. The position discrepancy, when the contact is made, is unavoidable whatever is the controller or the approach (unless a very prediction is made in the master side), this is due to the undergo physical time-delay. Nevertheless, the virtual probe position x e is stably maintained by the operator during the whole contact time. Several simulations are conducted with multiple hard and viscous contacts, they show that the behavior of the force feedback interaction is stable whatever the time delay. Obviously, one must not suspect that functional performances are acceptable for an actual use in the presence of important time delays.Figure5shows the result of the previously set up system under similar parameters and time-varying delay. The variable time delay τ 1 (t) is represented on the same figure, it satisfies the causality constraints: tτ 1 (t) > 0 and τ 1 (t) < 1 ∀t. In the simulation 0.5 6 τ 1 (t) 6 1.5sec. Three different simulations are

Figure 6 :

 6 Figure 6: Simulation behavior with an error estimation of the master model parameters (mass and friction).

  F e .I nt h i ss i m u l a t i o ni nfigure 6 m =0 .3 kg and b =5Nm/s are the estimated parameters, we can notice that the behavior of the interaction is still stable, figure 7.

Figure 7 :

 7 Figure 7: Stability margin expressing robustness of the control scheme

The developped proof holds for the non-linear haptic model described by the classical dynamic equation: M(q)q + C(q, q) q + G(q). The demonstration is also trivial, we choose the linear case for the clarity of the presentation.

[START_REF] Niemeyer | Stable adaptive teleoperation[END_REF] parameters α can be set to an adequate value to improve performances.


Sofia Zaidenberg

Franc ¸ois Brémond Bernard Boulay

A generic framework for video understanding applied to group behavior recognition

This paper presents an approach to detect and track groups of people in video-surveillance applications, and to automatically recognize their behavior. This method keeps track of individuals moving together by maintaining a spacial and temporal group coherence. First, people are individually detected and tracked. Second, their trajectories are analyzed over a temporal window and clustered using the Mean-Shift algorithm. A coherence value describes how well a set of people can be described as a group. Furthermore, we propose a formal event description language. The group events recognition approach is successfully validated on 4 camera views from 3 datasets: an airport, a subway, a shopping center corridor and an entrance hall.

Introduction

In the framework of a video understanding system (figure 1), video sequences are abstracted in physical objects: objects of interest for a given application. Then the physical objects are used to recognize events. In this paper, we are interested by the group behavior in public spaces. Given a set of detected and tracked people, our task is finding associations of those people into spatially and temporally coherent groups, and detecting events describing group behavior. Tracking people, and especially groups of people in relatively an unconstrained, cluttered environment is a challenging task for various reasons. In [START_REF] Ge | Automatically detecting the small group structure of a crowd[END_REF], Ge et al. propose a method to discover small groups of people in a crowd based on a bottom-up hierarchical clustering approach. Trajectories of pedestrians are clustered into groups based on their closeness in terms of distance and velocity. The experiments of this work have been made on videos taken from a very elevated viewpoint, providing few occlusions. Haritaoglu et al. [START_REF] Haritaoglu | Detection and tracking of shopping groups in stores[END_REF] detect and track groups of people as they shop in a store. Their method is based on searching strongly connected components in a graph created from trajectories of individual people, following the idea that people belonging to the same group have a lower inter-distance. This method however does not allow group members to move away and return to the group without being disconnected from it. Furthermore, the application of a shopping queue lacks genericity (people are rather static and have a structured behavior), it is not clear how well this method is adaptable to another context of use. Other approaches, such as [START_REF] Naturel | Detecting queues at vending machines: A statistical layered approach[END_REF], aim at detecting specific group-related events (e.g. queues at vending machines) without tracking. Here again, the method does not aim at consistently tracking a group as its dynamics vary. In [START_REF] Jacques | Understanding people motion in video sequences using voronoi diagrams[END_REF], an algorithm for group detection and classification as voluntary or involuntary (e.g. assembled randomly due to lack of space) is proposed. A top-down camera is used to track individuals, and Voronoi diagrams are used to quantify the sociological concept of personal space. No occlusion handling is done in this work hence the applicability to other points of view of the camera or to denser scenes is questionable. Figure 2 shows the result of our event recognition on tracked groups on a dataset recorded at the Eindhoven airport.

Event recognition is a key task in automatic understanding of video sequences. In this work we are mainly interested in group events, but the usual techniques can be applied to different kinds of objects (person, vehicle, group,...). The typical detection algorithm (figure 1) takes as input a video sequence and extracts interesting objects (physical objects). This abstraction stage is the layer be- tween the image and the semantic worlds. Then, these objects of interest are used to model events. Finally, the events are recognized. The abstraction stage determines which modeling techniques can be applied. The possible abstraction technique can be pixel based [START_REF] Bobick | The recognition of human movement using temporal templates[END_REF] or object based [START_REF] Vu | Automatic video interpretation: a novel algorithm for temporal scenario recognition[END_REF]. The first kind of techniques is not well adapted for groups. Indeed, persons belonging to the same group are not necessary physically linked. With object abstraction, a video sequence is represented thanks to the detected objects (persons, vehicles,...) and their associated properties (speed, trajectory,...). In the literature, several approaches use this abstraction level [START_REF] Romdhane | A framework dealing with Uncertainty for Complex Event Recognition[END_REF], because an activity can be naturally modeled based on those properties.

Lavee et al. [START_REF] Lavee | Understanding video events: A survey of methods for automatic interpretation of semantic occurrences in video[END_REF] classify existing event modeling techniques in three categories: the pattern recognition models, the state based models and the semantic models. The pattern recognition models are classical recognition techniques using classifiers as the nearest neighbor method, boost techniques, support vector machines and neural networks [START_REF] Chen | An interactive semantic video mining and retrieval platform-application in transportation surveillance video for incident detection[END_REF]. These techniques are well formalized. But adding new types of events implies the training of new classifiers. The state based models formalize the events in spatial and temporal terms using semantic knowledge: Finite State Machines (FSM), Bayesian Networks (BN), Hidden Markov Models (HMM), Dynamic Bayesian Networks (DBN) and Conditional Random Fields (CRF). The HMMs and all their variants are heavily used for event modeling [START_REF] Duong | Activity recognition and abnormality detection with the switching hidden semi-markov model[END_REF]. They take the advantages of the FSM (temporal modeling) and of the BNs (probabilistic modeling). But due to the nature of the HMMs (time sliced structure), the complex temporal relations (e.g. during) are not easily modeled. Lin et al. [START_REF] Lin | Group Event Detection With a Varying Number of Group Members for Video Surveillance[END_REF] propose an asynchronous HMM to recognize group events. Brdiczka et al. construct in [START_REF] Brdiczka | Automatic detection of interaction groups[END_REF] HMMs upon conversational hypotheses to model group events during a meeting. One drawback of the modified HMM methods is that since the classical structure of HMMs is modified, efficient algorithms can not be applied without approximation. The semantic models define spatio-temporal relations between sub-events to model complex events. Due to the nature of these models, the events must be defined by an expert of the application domain. Moreover, these models are often deterministic. Several techniques are studied: grammar based models, Petri nets (PN), constraint solving models and logic based models. As shown in this section, the quantity of techniques for abstraction and event modeling is huge. In this paper, we propose a framework (ScReK: Scenario Recognition based on Knowledge) to easily model the semantic knowledge of the considered application domains: the objects of interest and the scenario models, and to recognize events associated to the detected group based on spatio-temporal constraints.

In the rest of this paper, we first describe our technique to detect and track groups of people (section 2), then we describe our event detection method applied to tracked groups (section 3). Section 4 presents evaluations.

Group Tracking

Given a set of detected and tracked people, the proposed method focuses mainly on the task of finding associations of those people into spatially and temporally coherent groups. The human definition of a group is people that know each other or interact with each other. In fact, according to McPhail [START_REF] Mcphail | Using film to analyze pedestrian behavior[END_REF]: Two defining criteria of a group [are] proximity and/or conversation between two or more persons. It is quite difficult to directly detect people interactions and conversation in a video or the fact that people know each other. For automatic recognition we derive this definition: two or more people who are spatially and temporally close to each other and have similar direction and speed of movement, or better: people having similar trajectories.

Group tracking is based on people detection. The people detection can be performed by various methods. We have compared several methods and chosen the best one, it is based on background-subtraction described in [START_REF] Yao | Multi-layer background subtraction based on color and texture[END_REF] because of the quality of its results (see table 1 for a comparison of several methods).

Blobs of foreground pixels are grouped to form physical objects (also called mobiles) classified into predefined categories based on the 3D size of objects (using a calibrated camera): GROUP_OF_PERSONS, PERSON and NOISE. When people overlap (which happens quite often with a low viewpoint, such as in figure 6) or are too close to each other, segmentation fails to split them and they are detected as a single object classified as GROUP_OF_PERSONS because its size is bigger than the size of a single person. Thoses classes of objects are specified using gaussian functions. Mean, sigma, min and max values are provided for each class and a score is computed representing how well an object's dimensions fit in each category. The category with the best score is assigned as the class of the object. Detected objects at each frame are tracked consistently on the long term using a multiple feature-based tracker [START_REF] Chau | A multi-feature tracking algorithm enabling adaptation to context variations[END_REF].

Individual trajectories are the input of the group tracking algorithm, which is divided into four parts: creation, update, split/merge and termination. In order to detect temporally coherent groups, we observe people trajectories over a time window, denoted delay T . In the experiments presented section 4, we used T = 20 frames. Working at frame t c -T , t c being the current frame of the video stream, we cluster trajectories of individuals between frames t c -T and t c to find similar trajectories, representative of groups. We choose the Mean-Shift clustering algorithm [START_REF] Fukunaga | The estimation of the gradient of a density function, with applications in pattern recognition[END_REF] because it does not require to set as input the number of clusters. However, Mean-Shift does require a tolerance parameter determining the size of the neighborhood for creating clusters. Figure 3 shows the input prints and the clustering result.

A trajectory is defined as

T raj = {(x i , y i ), i = 0 . . . T -1}∪{(s xi , s yi ), i = 1 . . . T -1} where (x i , y i ), i ∈ [0; T -1]
in each trajectory is the position of a group in the same frame i, and

(s xi , s yi ) = speed(i-1, i), i ∈ [1; T -1]
is the speed of the group between frames i -1 and i. If k positions on the trajectory are missing because of lacking detections, we interpolate the k missing positions between known ones. Each trajectory is a point in a 2(2T -1)dimensional space. Mean-Shift is applied on a set of such points. To make the approach more generic and being able to add other features, we normalize the values using minimum and maximum ranges. The range of positions on the ground plane is determined by the field of view. The minimum speed is 0 and the maximum speed is set to 10 m/s, greatly exceeding all observed values. From the raw value of x, y and s (the speed) denoted by r ∈ [min, max], we compute the corresponding normalized value n ∈ [0, 1] as: n = r-min max-min , where min and max are the respective minimum and maximum values. We set the tolerance to 0.1, considering grouping trajectories distant by less than 10% of the maximum. This value is quite low because clustering is used only to group very close people, the case where people temporarily split being handled by the update step described below.

We characterize a group by three properties: the average over the frames in which the group is detected of the inter-mobile distance and the average over frames of standard deviations of speed and direction. These properties enable the definition of a coherence criterion: groupIncoherence = ω 1 • distanceAvg + ω 2 • speedStdDev + ω 3 • directionStdDev, where the weights ω 1 , ω 2 and ω 3 are normalization parameters. We use ω 1 = 7 and ω 2 = ω 3 = 5 to slightly favor distance over speed and direction similarity which are quite noisy. With this definition, a low value of groupIncoherence is significative of a group.

Groups are created from clusters of more than one physical object. In the case where one GROUP_OF_PERSONS object is detected at frame t c -T , we analyze its trajectory through the time window. If this object stays the size of a group, or is close to other objects, we can create a group and compute its groupIncoherence. If the resulting value is low enough, we keep the created group. In case of a single GROUP_OF_PERSONS object, the groupIncoherence value is naturally very low because of a null distanceAvg component. The creation step is made up of these two cases.

Group dynamics vary. Sometimes all group members do not have similar trajectories, for example when the group is waiting while one member buys a ticket at a vending machine. Clustering is not enough to correctly update an existing group in that case. First, we try to associate clusters with groups existing at the previous frame, using the notion of probable group of a mobile, defined hereafter. During tracking, mobiles detected at different frames are connected by probabilistic links in order to track consistently the same real objects. We use the term father and son for the mobiles resp. in the oldest and most recent frame of the link. If a father, within a window of T frames, of the mobile m was in a group g and the link probability between father and son is above a given threshold (a value of 0.6 is usually used in the experiments section 4), then the father's group g is called the probable group of the mobile m: P G(m) = g. Each cluster c is associated with the probable group of most mobiles in the cluster: G(c) = argmax g∈{g c i } |{g c i |g c i = g}|, where G(c) is the group associated with cluster c and {g c i } = {P G(m c i )} the set of probable groups of mobiles belonging to cluster c ({m c i } being the set of mobiles in cluster c). Several clusters can be associated to the same group, ensuring that group members having temporarily diverging trajectories will be kept in the group for a minimal amount of time. Each mobile m c i is added to the group G(c) if this group is really the probable group of the considered mobile: P G(m c i ) = G(c). In fact, the update step aims at tracking existing members of the group and not new comers. This procedure is summarized in algorithm 1.

The split of groups operates naturally. When a mobile from a group has moved away for too many frames, its probable group becomes empty and it cannot be added to an existing group during the update step, so it splits. It may be part of a new group in the creation step, if it gets clustered together with other mobiles.

Two groups g 1 and g 2 can be merged if two mobiles, one in each group at frame

t c -T + k (k ∈ [0; T -1]),
c i ∈ {m c } do g c i = P G(m c i ); G(c) = argmax g∈{g c i } |{g c i |g c i = g}|; for m c i ∈ mobiles tc-T do if P G(m c i ) = G(c) then G(c).add(m c i );
g 1 and g 2 is kept and all mobiles of the disappearing group are added into the remaining group. The group termination step erases old groups. Mobiles that have been detected at a largely outdated frame (e.g. t c -5T ) are deleted at frame t c -T and empty groups are erased. As a consequence, groups having no new mobiles for 5T frames are erased. All existing groups, even currently empty ones, can potentially be updated.

Finally, the output of the group tracker, which is the input of the event detection, is a set of tracked groups (keeping a consistent id through frames) having properties (such as the intra-objects distance) and composed of detected physical objects at each frame.

Event Recognition: a Generic Framework

In this work, a generic framework for event recognition is proposed (ScReK). The motivation of this framework is to use it within any video understanding application. The genericity is obtained in terms of objects of interest and event models. We can identify two main parts in an event recognition process: the application knowledge (what are the expected objects? what are the event models?) and the event recognition algorithm. Knowledge representation is a key issue for genericity. We believe that the knowledge should be modeled with the collaboration of two different categories of people (figure 4): vision experts (specialists in vision algorithms), and application domain experts (specialists in the expected events of their domain). Vision experts are modeling the objects of interest (detected by vision algorithms) and the video primitives (properties computed on the detected objects of interest). Domain experts have to model the expected application events. Usually, for video event recognition, knowledge is represented using OWL (Web Ontology Language). Even with tools like Protégé, it is difficult for a non computer specialist to create her/his own model without a long and tedious learning of the OWL formalism. The ScReK framework proposes its own declarative language to easily describe the application domain knowledge: the ontology. ScReK proposes a grammar description of the objects and events using the extended BNF (Backus Naur Form) representation. O i , is described by its parent, O j , and its attributes:

O i = {a k } k=0,...O n i .
The objects are defined using an inheritance mechanism. The object O i inherits all the attributes of its parent O j . The attributes are described with the help of basic types. 11 basic types are predefined: boolean, integer, double, timestamp, time interval, 2D point (integer and double), 3D point (integer and double), and list of 3D points. The user can contribute by adding new basic types. Moreover, a history of the values of a basic type is automatically stored. It is useful for vision primitives based on the evolution of a value in time (e.g. trajectory). For group behavior recognition, detected group objects within the video sequence and scene context objects (zone, equipment) are described. The scene context objects help to recognize specific events (e.g. by defining a forbidden access zone or a threshold). For instance, the class of group objects is defined as follows in the ScReK language: A Group is a Mobile and it inherits all the attributes of a Mobile object (3D size, 3D position,...). A Group is not constant (dynamic, i.e. its attributes values can change throughout time). One of its attributes, NumberOfMobiles is the number of objects which compose the group. The second kind of knowledge to represent is the event models. They are composed of 6 parts: (1) the type of the scenario can be one of the following: PrimitiveState, Com-positeState, PrimitiveEvent, CompositeEvent, from the simplest to the most complex events. ( 2) the name of the event model which can be referenced for more complex events.

(3) the list of physical objects (i.e. objects of interest) involved in the event. The type of the objects is depending on the application domain. (4) the list of components contains the sub-events composing the event model. ( 5) the list of constraints for the physical objects or the components. The constraints can be temporal (between the components) or symbolic (for physical objects). ( 6) the alarm information describes the importance of the scenario model in terms of urgency. Three values are possible, from less urgent to more urgent: NOTURGENT, URGENT, VERYURGENT.

The alarm level can be used to filter recognized events, for displaying only important events to the user. Hereafter is a sample event model: C o m p o s i t e E v e n t ( b r o w s i n g , P h y s i c a l O b j e c t s ( ( g : Group ) , ( e : E q u i p m e n t ) ) Components ( ( c1 : G r o u p S t o p ( g ) ) ( c2 : G r o u p N e a r E q u i p m e n t ( g , e ) ) ) C o n s t r a i n t s ( ( e->Name = " shop window " ) ) Alarm ( ( L e v e l : URGENT ) ) )

The application domain expert models the event browsing by "a group is stopped in front of the shop-window" with the model above. The vision expert models the subevent Group Near Equipment (by measuring the distance between a group and an equipment) and Group Stop (by computing the speed of a group).

The last part of the event recognition framework is the recognition algorithm itself. The proposed algorithm solves spatio-temporal constraints on the detected groups. The usual algorithms to recognize such events can be time consuming. The ScReK framework proposes to define optimal event models: at most two components, at most one temporal constraint (Allen's algebra) between these components. This property is not restrictive since all event models can be optimized in this format. Thanks to the optimal property, the event model tree is computed. The tree defines which sub-event (component) triggers the recognition of which event: the sub-event which happens last in time triggers the recognition of the global event. For instance, the event A has two components B and C with constraint: B before C. The recognition of C triggers the recognition of A. The tree triggers the recognition of the only events that can happen, decreasing the computation time. The first step of the event recognition process is to recognize all the possible simple events (most of these events are based on the vision primitives) by instantiating all the models with the detected objects (e.g. instantiating the model Group_Stays_Inside_Zone (takes as input one group and one zone) for all the detected groups and all the zones of the context). The second step consists in recognizing complex events according to the event model tree and the simple events previously recognized. The final step checks if the recognized event at time t has been already recognized previously to update the event (end time) or create a new one.

Results

People detection is an input to group detection. We compared several methods to validate our choice of method. Table 1 sums up the results of an evaluation done on a 36006 frames sequence (approximately 2 hours of video) in which 37 ground truth (GT) objects (people) have been annotated. [START_REF] Corvée | Haar like and LBP based features for face, head and people detection in video sequences[END_REF] is a feature-based people detector whereas [START_REF] Nghiem | Controlling Background Subtraction Algorithms for Robust Object Detection[END_REF] and [START_REF] Yao | Multi-layer background subtraction based on color and texture[END_REF] Figure 5. Proposed group event ontology. [START_REF] Corvée | Haar like and LBP based features for face, head and people detection in video sequences[END_REF] [ both perform motion segmentation and classification of detected objects. The method C combines the first two methods for a more robust detection than each one separately. The method from [START_REF] Yao | Multi-layer background subtraction based on color and texture[END_REF] gives the best results and is used as input of the group tracking process. This method learns a background model, resulting in better motion segmentation and better detection of small objects (far from the camera) and static objects. The drawback is the time necessary to learn the model and the low speed of the backgroundsubtraction.

We have performed evaluation of the group tracking algorithm using 4 different views from 3 datasets: videos recorded for the european project VANAHEIM1 in the Turin subway (figure 6), videos recorded for the european project ViCoMo2 at the Eindhoven airport (figure 2) and videos from the benchmarking CAVIAR3 dataset: the INRIA entrance and the shopping center corridor. In tables 2 and 3 the following metrics are used. The fragmentation metric computes throughout time how many tracked objects are associated with one reference object (ground truth data). The tracking time metric measures the percentage of time during which a reference object is correctly tracked. The purity computes the number of reference object IDs per tracked object. A value close to 1 is significative of a good tracking.

Table 2 shows group detection and tracking results on 16 sequences from the CAVIAR dataset. The first 9 are sequences from INRIA, and the remaining are from the shopping center corridor. One can notice that the shopping view is far more challenging than the hall because more people are visible and there are more occlusions due to the low position of the camera. Table 3 contains the results of this evaluation on 3 annotated sequences (resp. 128, 1373 and 17992 frames) from the Turin subway dataset. In both tables, detection results are good for almost all sequences. In the sequence c2ls1, ground truth groups in the far end of the corridor fail to be detected because of the limitations of the background-subtraction method. Tracking shows good results with a few exceptions. For instance, sequence 2 of table 3 contains a main group present in the foreground for the whole duration of the sequence. This group is correctly tracked with only one id-switch, but many groups are annotated far in the background and are difficult to detect for the motion segmentation algorithm. Their sparse detection results in many id-switches for group tracking. At the best of our knowledge, there is no possibility of comparing our method to an existing one (no public results or code available).

One major achievement of this paper is an ontology for group events based on video sensor (figure 5). The ontology is composed of 49 event models (45 models are generic and re-usable in any application with groups (Group stop, Group lively,...), 4 models are specifically defined for the applications of this paper (the events depend on the application context, enter shop,...)). The events have been modeled with help of metro surveillance staff.

Conclusions

We propose a generic, plug and play framework for event recognition from videos: ScReK. The scientific community can share a common ontology composed of event models and vision primitives. We demonstrate this framework on 4 group behavior recognition applications, using a novel group tracking approach. This approach gives satisfying results even on very challenging datasets (numerous occlusions and long duration sequences) such as in figure 6. The vision primitives are based on global attributes of groups (position, speed, size). The proposed event detection approach correctly recognizes events but shows its limitation for some specific events (e.g. fighting is best characterized by internal group movement). Adapted vision primitives, such as optical flow, solve specific limitations and are easy to plug into ScReK. Moreover, in this work the gap between video data and semantical events is modeled manually by vision experts, the next step is to learn automatically the vision primitives.

Figure 1 .

 1 Figure 1. Description of the proposed video understanding system

Figure 2 .

 2 Figure 2. Group event recognition on Eindhoven airport sequences. Left: The group is detected as splitting into 2 sub-groups. Right: Two groups are detected as merging into one.

Figure 3 .

 3 Figure 3. Example of 5 tracked people clustered into 3 trajectory clusters. A group is created from cluster 2.

1 :

 1 have the same son at frame t c -T + l (l ∈ [k + 1; T -1]), meaning that the two mobiles will merge. The oldest group among Algorithm Update of groups. input : {groups tc-T -1 }, {mobiles tc-T } output: updated {groups tc-T } {clusters tc-T } = M eanShif t({mobiles tc-T }); for c ∈ {clusters tc-T } do for m

Figure 4 .

 4 Figure 4. Knowledge modeling for video event recognition.

  c l a s s Group : M o b i l e { c o n s t f a l s e ; C S I n t NumberOfMobiles ; CSDouble A v e r a g e D i s t M o b i l e s ; }

Figure 6 .

 6 Figure 6. Event detection in Turin subway.

  instances of each event because we only focus on meaningful group events. The events are correctly recognized with low false positive and false negative rates. Most of the false positive detections for the event getting off train are due to the fact that the door in the foreground is detected as a person when open. The errors can be corrected by adding a new video primitive: door detector.

Figure 7 .

 7 Figure 7. Group event recognition for the CAVIAR sequences. a. Fighting then splitting. b. Exit from the shop. c. Browsing. d. The mis-detected group (ghost due to reflections) is browsing.

Table 1 .

 1 Comparison several people detection methods.

	15]	C	[18]

Table 2 .

 2 The results of the group event recognition are given in table4for the interesting events. Examples of event recognition are shown in figures 2, 6 and 7. There is only a few Results of group detection and tracking on 16 CAVIAR sequences. (Seq -official sequence name, Prec -Precision, Sens -Sensitivity, Frag -Fragmentation, TT -Tracking Time)

				Detection			Tracking
	Seq	TP	FP	FN	Prec	Sens	Frag	TT	Purity
	fc	125 101 3	0.55 0.98 1	0.97 1
	fomd 159 0	8	1	0.95 1	0.95 1
	fra1 139 0	67	1	0.67 1	0.61 1
	fra2 141 0	55	1	0.72 1	0.70 1
	mc1 231 0	82	1	0.74 1	0.72 1
	ms3g 145 37	37	0.80 0.80 1	0.61 1
	mwt1 156 0	89	1	0.64 1	0.28 1
	mwt2 336 0	268 1	0.56 1	0.53 1
	sp	165 4	36	1	0.82 1	0.67 1
	c2es1 858 652 487 0.57 0.64 0.58 0.41 0.81
	c2es3 1093 550 735 0.66 0.60 0.66 0.38 0.83
	c2ls1 788 1664 655 0.32 0.55 0.50 0.13 1
	c3ps1 1298 1135 210 0.54 0.86 1	0.68 1
	cosme2 1119 852 35	0.57 0.97 0.25 0.60 1
	csa1 269 163 0	0.63 1	1	0.96 1
	cwbs12224 89	1090 0.96 0.67 1	0.45 0.80
				Detection			Tracking
	Seq	TP	FP	FN	Prec	Sens	Frag	TT	Purity
	1	65	0	6	1	0.91 1	0.86 1
	2	1346 69	318 0.95 0.80 0.80 0.14 0.91
	3	6977 1677 4594 0.80 0.60 0.40 0.32 0.98

Table 3 .

 3 Results of group detection and tracking on 3 sequences from the Turin subway. (Prec -Precision, Sens -Sensitivity, Frag -Fragmentation, TT -Tracking Time)

		GT TP FP FN
	fighting	2	1	0	1
	split up	3	3	0	0
	joining	3	3	0	0
	shop enter	5	5	0	0
	shop exit	6	6	1	0
	browsing	3	3	1	0
	getting off train 10	9	8	1

Table 4 .

 4 Group event recognition for the 3 video datasets

https://www.vanaheim-project.eu/

http://www.vicomo.org/

http://homepages.inf.ed.ac.uk/rbf/CAVIARDATA1/

Acknowledgment This work was supported partly by the Video-Id, ViCoMo, Vanaheim, and Support projects. However, the views and opinions expressed herein do not necessarily reflect those of the financing institutions.