B Staber

J Guilleminot Ůniversité

Stochastic hyperelastic constitutive laws and identification procedure for soft biological tissues with intrinsic variability

Keywords: biological material, calibration, constitutive law, hyperelasticity, information theory, variability 2010 MSC: 74A40, 74B20, 74Q15, 74S60 1. Preliminaries

Stochastic hyperelastic constitutive laws and identification procedure for soft biological tissues with intrinsic variability

B Staber, J Guilleminot

Introduction

The constitutive modeling and data-driven calibration of the mechanical behavior of soft biological tissues is of paramount importance in a wide range of applications from bioengineering and computational -and potentially, multiscale -mechanics of hierarchical biomaterials (see [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF][START_REF] Holzapfel | Constitutive modelling 455 of arteries[END_REF][START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF] for the modeling of arteries, for instance) to computed-assisted surgery (CAS). While the definition of suitable models is still the matter of active research, due to the great complexity and multiphysics aspects of the mechanisms governing the macroscopic behavior, one additional difficulty lies in the fact that the experimental characterization of such materials typically exhibits a large variability; see, in a non-exhaustive manner, [START_REF] Miller | Constitutive modelling of brain tissue: experiment and theory[END_REF][START_REF]Mechanical properties of brain tissue in tension[END_REF][START_REF] Velardi | Anisotropic constitutive equations and experimental tensile behavior of brain tissue[END_REF][START_REF] Hrapko | Characterisation of the mechanical behaviour of brain tissue in 460 compression and shear[END_REF][START_REF] Rashid | Mechani-500 cal characterization of brain tissue in compression at dynamic strain rates[END_REF][START_REF] Rashid | Experimental 510 characterisation of neural tissue at collision speeds[END_REF][START_REF] Jin | A comprehensive experimental study on material properties of human brain tissue[END_REF][START_REF]Mechanical characterization of brain tissue in simple shear at dynamic strain rates[END_REF][START_REF]Mechanical characterization of brain tissue in tension at dynamic strain rates[END_REF], [START_REF] Sakuma | In vitro measurement of mechanical properties of liver tissue under compression and elongation using a new test piece holding method with surgical glue[END_REF]6,[START_REF] Roan | The nonlinear material proper-515 ties of liver tissue determined from no-slip uniaxial compression experiments[END_REF][START_REF]Strain rate-dependent viscohyperelastic constitutive modeling of bovine liver tissue[END_REF], [START_REF] Sparrey | Compression behavior of porcine spinal cord white matter[END_REF] and [START_REF] Umale | Experimental mechanical characterization of abdominal organs: liver, kidney & spleen[END_REF] for experimental characterizations of brain tissues, liver tissues, spinal cord white matter and abdominal organs, respectively. These uncertainties, which may be due to various ˚Corresponding author Email address: johann.guilleminot@u-pem.fr (J. Guilleminot) factors (such as microstructural randomness, subject-to-20 subject variability, age, gender or even potential disease), may critically impact simulation predictiveness or clinical diagnosis for example. In practice, these fluctuations are usually ignored and smoothed out through an averaging procedure along the 25 statistical dimension. This mean response is subsequently fitted by solving a least-squares problem by means of a deterministic optimization method (see [START_REF] Ogden | Fitting hyperelastic models to experimental data[END_REF]). Whereas such approximations may be found satisfactory in some cases, depending on the application under consideration, 30 performing high-fidelity computations or making informed decisions requires the modeling and calibration issues to be both addressed in a probabilistic framework. From a stochastic standpoint, the Bayesian paradigm provides a sound approach to data-driven modeling for deterministic [START_REF] Sakuma | In vitro measurement of mechanical properties of liver tissue under compression and elongation using a new test piece holding method with surgical glue[END_REF] and random parameters [START_REF] Bernado | Bayesian Theory[END_REF][START_REF] Gelman | Bayesian 445 Data Analysis[END_REF]. This natural path was very recently pursued in [START_REF]Bayesian calibration of hyperelastic constitutive models of soft tissue[END_REF] and previously used by the same authors for model selection in [START_REF] Madireddy | A bayesian ap-475 proach to selecting hyperelastic constitutive models of soft tissue[END_REF].

In this work, we propose an alternative methodology where the probability density functions defining the set 40 of randomized parameters are constructed explicitly. The main advantage of such a formulation (in contrast to a Bayesian approach where the posterior distributions are generally unknown analytically) is that it allows existence results for the solution of the associated boundary value problem to be invoked (with probability one). In practice, the proposed class of stochastic stored energy functions can purposely be used as:

• an admissible model making the identification task well-posed (from a mechanical standpoint) whenever experimental data are limited;

• a prior model in a Bayesian setting, especially when the number of the parameters defining the model (which are referred to as hyperparameters) is large; see [START_REF] Soize | A computational inverse method for identification of non-gaussian random fields using the bayesian approach in very high dimension[END_REF] for an example in three-dimensional linear elasticity.

The aim of this research is twofold. First, it attempts to lay out a framework, based on information theory, for the modeling of stochastic stored energy functions. To this purpose, we provide a self-contained and applicationoriented presentation of a stochastic model derived for the Ogden class of stored energy functions [START_REF] Staber | Stochastic modeling of a class of stored energy functions for incompressible hyperelastic materials with uncertainties[END_REF][START_REF]Stochastic modeling of the ogden class of stored energy functions for hyperelastic materials: the compressible case[END_REF]. This class turns out to be relevant to the modeling of some isotropic biological materials (see [START_REF] Chagnon | Hyperelastic energy densities for soft biological tissues: A review[END_REF] for example). Additionally, we address the definition of a calibration methodology. The latter specifically allows for the identification of the model hyperparameters using mean and variance estimates along a loading curve. Second, it aims at assessing the capability of these models to reproduce the experimental variability that is observed in the mechanical response of biological tissues at large strains (and small strain rates). From a clinical perspective, these statistical fluctuations may be considered as informative quantities in the early stage prognosis of diseases affecting the mechanical response of the tissues under consideration.

This outline of this paper is as follows. In Section 2, we first set up notation and provide a brief reminder on deterministic constitutive modeling in hyperelasticity. The information-theoretic stochastic model and the calibration procedure are then presented in Section 3. Finally, some applications involving experimental results taken from the literature are discussed in Section 4.

Notation

Throughout this paper, deterministic scalars, vectors and second-order tensors are denoted by a, a and rAs. which may change from line to line.

Their stochastic counterparts are respectively denoted by

Deterministic constitutive model

Basic notation and kinematics

Let Ω be the region in R 3 occupied by the body under consideration at rest and let BΩ denote its boundary. The 95 body is assumed to undergo a deformation ϕ : Ω Q x Þ Ñ x ϕ P Ω ϕ , where Ω ϕ is the current configuration. The mapping ϕ is classically assumed to be sufficiently smooth, orientation-preserving and injective.

The associated deformation gradient rF s is defined as

rF pxqs " ∇ x ϕpxq , @x P Ω , (1)
and satisfies the condition detprF sq ą 0 in Ω. For the sake of notational convenience, spatial indexing is omitted from now on. The principal stretches of the deformation gradient are denoted as tυ j prF squ

where rT s is the first Piola-Kirchhoff stress tensor (which is also referred to as the engineering stress), rF s denotes the deformation gradient (with unitary determinant) and is a Lagrange multiplier raised by the incompressibility condition. Additionally, we focus on Ogden-type stored energy functions for isotropic hyperelastic materials. Such functions are defined, at orders m ě 1 and n ě 1, as [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF] wprF sq "

m ÿ i"1 p i Φ ηi prF sq `m`n ÿ j"m`1 p j Υ ηj prF sq , (4) with
Φ ηi prF sq " υ ηi 1 `υηi 2 `υηi 3 ´3 (5)
(with a slight abuse of notation) and Υ ηj prF sq " pυ 1 υ 2 q ηj `pυ 1 υ 3 q ηj `pυ 2 υ 3 q ηj ´3 ,

where it is recalled that tυ j u 3 j"1 are the principal stretches associated with the deformation gradient rF s.

Upon expanding the stored energy function near the reference configuration, it is found that the parameters p " pp 1 , . . . , p np q and η " pη 1 , . . . , η np q must satisfy the following additional consistency relation [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF]:

np ÿ k"1 p k η 2 k " 2µ , (7)
where µ ą 0 is the shear modulus.

The applications considered in this paper consist in simple compression tests along a principal direction, say e 3 " p0, 0, 1q. In this case, and by combining Eq. ( 3) and the incompressibility condition υ 1 υ 2 υ 3 " 1, one can show (see e.g. [START_REF] Ogden | Non-Linear Elastic Deformations[END_REF][START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF]) that the resulting principal Piola-Kirchhoff stress σ writes σpυq "

m ÿ i"1 p i η i pυ ηi´1 ´υ´1´ηi{2 q ´m`n ÿ j"m`1 p j η j pυ ´ηj ´1 ´υ´1`ηj{2 q , ( 8
)
where the substitution υ " υ 3 has been introduced for notational convenience.

Fundamental assumptions

From a theoretical standpoint, the available existence theorems for boundary value problems in hyperelasticity (see the references hereinafter) rely on the existence of minimizers for the total potential energy function and were 110 established by resorting on the calculus of variations. Not surprisingly, the investigation of sufficient conditions on the stored energy function such that the aforementioned existence results can effectively be proven has attracted much attention from researchers in applied mathematics and theoretical mechanics. To date, minimal requirements are essentially concerned with some generalized convexity conditions, as well as with a coerciveness inequality that models an asymptotic behavior at large strains; see the seminal work reported in [START_REF] Morrey | Quasi-convexity and the lower semicontinuity of multiple integrals[END_REF], as well as the textbook [START_REF] Ciarlet | Mathematical elasticity: Three-dimensional elasticity[END_REF].

Certainly one of the most invoked assumption is that of polyconvexity, introduced by Ball in [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF]. This fundamental property, when complemented by a set of suitable growth conditions, insures the existence of at least one minimizer for the total potential energy; hence ensuring the well-125 posedness of the nonlinear boundary value problem (see also [START_REF] Balzani | A polyconvex framework for soft biological tissues. adjustment 425 to experimental data[END_REF] for a specific discussion on biological materials).

In what follows, and following the result by Müller et al. [START_REF] Müller | On a new class of elastic deformations not allowing for cavitation[END_REF], we assume that (see [START_REF]Some recent developments in nonlinear elasticity and 420 its applications to materials science[END_REF])

p k ą 0 , 1 ď k ď n p , (9a)
η 1 ě ¨¨¨ě η m ě 1 , η 1 ě 2 , (9b)
η m`1 ě ¨¨¨ě η np ě 1 , η m`1 ě 3 2 , (9c)
with n p " m `n, ensuring the existence of at least one minimizer for the total potential energy.

Remark 1. The single term stored energy function

wprF sq " p pυ η 1 `υη 2 `υη 3 ´3q (10)
is covered by the existence results if p ą 0 and η ě 3 (see [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF]). Moreover, the consistency with linearized elasticity yields pη 2 " 2µ, so that the stored energy function writes:

wprF sq " 2µ η 2 pυ η 1 `υη 2 `υη 3 ´3q . ( 11
)
The particular case η " 2 leads to the well known Neo-130

Hookean model in three-dimensional elasticity, which is not covered by the existence results invoked in this paper.

In the next section, we first discuss the randomization of Ogden-type stored energy functions. A new methodology allowing for the calibration of the probabilistic model is 135 subsequently proposed.

Stochastic hyperelastic constitutive models

The generalization of hyperelastic constitutive laws in a stochastic context has received little attention to date. The main difficulty raised by such a generalization stems 140 from the fact that the stochastic stored energy functions are simultaneously required (i) to accurately reproduce the experimental trends, and (ii) to be consistent with the theory of nonlinear elasticity (see Section 2.3). It should be noticed that considering one single constraint apart from 145 the other one may lead to models that • fit dispersed data very well (in terms of mean and variance at a given stretch for example), but are not consistent with theoretical aspects -hence yielding non-admissible samples for subsequent analysis; 150

• are fully consistent with existence theorems, but may not satisfactorily reproduce observed experimental results.

Clearly, any of the above situations reduces the predictive capabilities of the associated modeling framework. Note 155 in addition that experimental data are often very limited, so that the direct computation of statistical estimators at a reasonable level of convergence is itself an issue -in this case, the use of the maximum likelihood method may alternatively be considered.

160

Very recently, new contributions accommodating the aforementioned modeling constraints in the framework of information theory [START_REF] Jaynes | Information theory and stastitical mechanics ii[END_REF][START_REF]Information theory and statistical mechanics i[END_REF] were proposed; see [START_REF] Staber | Stochastic modeling of a class of stored energy functions for incompressible hyperelastic materials with uncertainties[END_REF] for the case of incompressible materials, as well as [START_REF]Stochastic modeling of the ogden class of stored energy functions for hyperelastic materials: the compressible case[END_REF] for the compressible case. In these works, stochastic stored energy 165 functions are constructed by randomizing, in a proper way, some parameters involved in the definition of classical, e.g. Ogden-type stored energy functions. More specifically, and following the maximum entropy (MaxEnt) principle, the probability density functions defining these random 170 parameters are constructed by maximizing the so-called Shannon differential entropy [START_REF] Shannon | A mathematical theory of communication[END_REF] over the set of probability density functions that satisfy a set of contraint equations (see Appendix A for more details). The latter, which are expressed in the form of mathematical expectations, are essentially induced by the information that is available on the set of random model parameters. This information is typically inferred from the properties ensuring the wellposedness of the associated variational formulation, such as positiveness, or from known statistical moments. By construction, such an approach reduces the modeling bias, especially when compared to methodologies relying on the a priori choice of standard statistical distributions, and readily ensures that all sufficient mathematical conditions are fulfilled (in accordance with the results from Section 2.3). The main result derived in [START_REF] Staber | Stochastic modeling of a class of stored energy functions for incompressible hyperelastic materials with uncertainties[END_REF] is first presented for completeness in Section 3.1. A novel methodology allowing for the calibration of the stochastic stored energy function using typically available experimental results (namely, the mean and variance evolutions) is then proposed in Section 3.2. This calibration strategy will be tested for various biological tissues in Section 4.

Construction of stochastic Ogden-type stored energy functions

Let P and µ be the random variables corresponding to the stochastic modeling of the vector-valued parameter p and the shear modulus µ, respectively. The stochastic version W of the Ogden potential w is then written as

W prF sq " m ÿ i"1 P i Φ ηi prF sq `m`n ÿ j"m`1 P j Υ ηj prF sq , ( 12
)
where Φ ηi and Υ ηj are defined by Eqs. ( 5)-( 6), and the exponents tη k u np k"1 are kept deterministic and such that the restrictions given by Eqs. (9b)-(9c) hold. The random variables P and µ also satisfy the consistency relation (see Eq. ( 7))

np ÿ k"1 P k η 2 k " 2µ (13)
almost surely. It should be noticed that the stochastic models of P and µ completely define the probabilistic model of the stochastic stored energy function W for a given value of the deformation gradient rF s.

In a first step, the probability density function p µ of the random shear modulus µ is constructed by taking into account that µ ą 0 almost surely, the mean constraint

E tµu " µ , µ ą 0 , (14)
and the additional constraint:

E tlogpµqu " ν , |ν| ă `8 . ( 15
)
From a theoretical standpoint, Eq. ( 15) ensures that µ and µ ´1 are both second-order random variables (see [START_REF] Guilleminot | On the statistical dependence for the components of random elasticity tensors exhibiting material symmetry properties[END_REF]).

According to the MaxEnt principle (see Appendix A.2 for technical details), the random variable µ then follows a Gamma distribution:

µ " Gpρ 1 , ρ 2 q , (16)
where the parameters ρ 1 ą 0 and ρ 2 ą 0 are raised by the MaxEnt principle.

200

In a second step, we introduce for technical convenience an auxiliary random variable U " pU 1 , . . . , U np´1 q such that (see Eq. ( 13))

P k " 2µ η 2 k U k , 1 ď k ď n p ´1 (17)
and

P np " 2µ η 2 np ˜1 ´np´1 ÿ k"1 U k ¸. ( 18
)
By construction (see Eq. (9a)), the random variable U takes its values in the admissible set

S " tu P ps0, 1rq np´1 | 1 ´řnp´1 k"1 u k ą 0u ( 19
)
that corresponds to the interior of a pn p ´1q-dimensional simplex. The probability density function p U of U is derived by invoking the principle of maximum entropy under constraints that are induced by Eq. (9a) (or equivalently, by Eq. [START_REF] Luenberger | Optimization by vector space methods[END_REF]; see [START_REF] Staber | Stochastic modeling of a class of stored energy functions for incompressible hyperelastic materials with uncertainties[END_REF]), namely:

E tlog pU k qu " ν k , 1 ď k ď n p ´1 , (20)
E # log ˜1 ´np´1 ÿ k"1 U k ¸+ " ν np , (21)
where |ν k | ă `8 for 1 ď k ď n p . It can then be shown (see Appendix A.3) that the random variable U follows a Dirichlet distribution:

U " Dirpλ 1 , . . . , λ np q , (22)
where the strictly positive parameters pλ 1 , . . . , λ np q are the Lagrange multipliers that are introduced by the MaxEnt principle.

Note that the marginal probability distributions of U thus obtained are unimodal when λ k ą 1, 1 ď k ď n p .

205

Whenever required, the case of multimodal distributions can be handled, either by defining a mixture of the proposed unimodal distributions (provided that a one-to-one correspondence between the modes of the random variables can be identified), or by using other probabilistic re-210 presentations, such as mixtures of polynomial chaos expansions (PCEs) [START_REF] Nouy | Identification of multi-modal random variables through mixtures of polynomial chaos expansions[END_REF] or PCEs defined with respect to multimodal probability distributions [START_REF]Polynomial chaos expansion of a multimodal random vector[END_REF].

The main result derived in [START_REF] Staber | Stochastic modeling of a class of stored energy functions for incompressible hyperelastic materials with uncertainties[END_REF] then states that the stochastic stored energy function defined by Eq. ( 12), where

215

• P is defined by Eqs. ( 17)-( 18),

• µ " Gpρ 1 , ρ 2 q, with ρ 1 ą 0 and ρ 2 ą 0,

• U " Dirpλ 1 , . . . , λ np q with λ i ą 0 for 1 ď i ď n p , is covered by the existence results (outlined in Section 2.3) and consistent with linearized elasticity almost surely. It is therefore admissible almost surely, which in turn implies that the associated stochastic nonlinear boundary value problem admits at least one solution.

The stochastic constitutive law is finally defined as rT s " BW prF sq BrF s

´ rF s ´T ,

where rT s denotes the random first Piola-Kirchhoff stress tensor. Using Eq. ( 8), it can be deduced that the random Piola-Kirchhoff stress Σ for a simple tensile (or compression) test writes Σpυq "

m ÿ i"1 P i η i pυ ηi´1 ´υ´1´ηi{2 q ´m`n ÿ j"m`1 P j η j pυ ´ηj ´1 ´υ´1`ηj{2 q . ( 24
)
Remark 2. It is observed that the variance of the Piola-Kirchhoff stress Σpυq is monotonically increasing with υ.

It follows that the proposed stochastic approach is indeed relevant to situations where the variance increases along the loading path (this feature being frequently encountered in practice). Nevertheless, some cases where the variance of the experimental stress does not increase together with the principal stretch may be addressed by randomizing the model exponents tη k u np k"1 . While such a model could certainly be constructed by pursuing a similar informationtheoretic methodology, such a strategy would introduce a set of hyperparameters belonging to high-dimensional spaces -hence making the calibration step more intricate. Remark 3. A robust generator for the random variable U can be obtained by introducing a set of Gamma random variables tY i u np i"1 , with

Y i " Gpλ i , 1q , (25)
and by sampling U through the component-wise mapping (see [START_REF] Devroye | Non-Uniform Random Variate Generation[END_REF])

U i :" Y i ˆ˜np ÿ k"1 Y k ¸´1 , 1 ď i ď n p ´1 . (26)
In the case n p " 2, U becomes a scalar random variable U which follows a Beta distribution. The Matlab functions gamrnd and betarnd might be useful in order to generate samples of Gamma and Beta random variables.

Identification of the stochastic stored energy functions based on experimental data: methodological aspects

In this section, we propose a methodology aiming at the calibration of the model parameters, which are gathered in the vector-valued Lagrange multipliers λ " pλ 1 , . . . , λ np q 245 and ρ " pρ 1 , ρ 2 q, using standard experimental data.

More specifically, let tξ exp k u N k"1 and tς exp k u N k"1 denote the sets of values for the mean and standard deviation of the random stress at prescribed stretches tυ exp k u N k"1 . These sets define the experimental information based on which 250 the stochastic representation has to be identified. For ease of explanation, we consider below the case of first-order stored energy functions (m " n " 1; hence n p " 2). The approach can easily be generalized to arbitrary orders at the expense of notational complexity. Closed-form expres-255 sions that are relevant to this case are first derived in Section 3.2.1. The two sequential steps involved in the calibration methodology are then described in Sections 3.2.2 and 3.2.3.

3.2.1.

Closed-form relations for first-order stochastic stored 260 energy functions When m " n " 1, the random stored energy function writes

W prF sq " P 1 pυ η1 1 `υη1 2 `υη1 3 ´3q `P2 ppυ 1 υ 2 q η2 `pυ 1 υ 3 q η2 `pυ 2 υ 3 q η2 ´3q , (27)
where the random variables P 1 and P 2 are given by

P 1 " 2µ η 2 1 U , P 2 " 2µ η 2 2 p1 ´U q , (28)
with η 1 ě 2 and η 2 ě 3{2. It is recalled that µ follows a Gamma distribution with parameters pρ 1 , ρ 2 q, and that the random variable U follows a Dirichlet distribution with parameters pλ 1 , λ 2 q. According to Eq. ( 24), the associated random Piola-Kirchhoff stress Σ is given, for a uniaxial tensile (or compression) test, by

Σpυq "P 1 Ψ 1 pυq `P2 Ψ 2 pυq , (29)
where the functions Ψ 1 and Ψ 2 are defined as

Ψ 1 pυq " η 1 ´υη1´1 ´υ´1´η1{2 ¯(30) and Ψ 2 pυq " ´η2 ´υ´η2´1 ´υ´1`η2{2 ¯. (31)
Moreover, it is readily deduced from Eq. ( 29) that the variance of Σpυq is given by

VartΣpυqu "VartP 1 uΨ 2 1 pυq `VartP 2 uΨ 2 2 pυq `2CovtP 1 , P 2 uΨ 1 pυqΨ 2 pυq , (32)
where VartP j u and CovtP 1 , P 2 u denote the variance of P j and the covariance between P 1 and P 2 , respectively. Next, and upon using the probabilistic properties of µ and U , it is found that:

VartP 1 u " 4ρ 1 ρ 2 2 η 4 1 λ 1 pλ 1 `λ2 `λ1 λ 2 `λ2 1 `ρ1 λ 2 q pλ 1 `λ2 q 2 p1 `λ1 `λ2 q , (33)
VartP 2 u "

4ρ 1 ρ 2 2 η 4 2 λ 2 pλ 1 `λ2 `λ1 λ 2 `λ2 2 `ρ1 λ 1 q pλ 1 `λ2 q 2 p1 `λ1 `λ2 q , (34)
and

CovtP 1 , P 2 u " 4ρ 1 ρ 2 2 η 2 1 η 2 2 λ 1 λ 2 pλ 1 `λ2 ´ρ1 q pλ 1 `λ2 q 2 p1 `λ1 `λ2 q . ( 35
)
Substituting Eqs. ( 33), ( 34) and ( 35) into Eq. ( 32) allows us to estimate the variance of the stress Σpυq for given values of the hyperparameters (without resorting to Monte Carlo simulations).

Defining Step 1: calibration of the mean model

In the first step of the approach, which is formulated deterministically, the mean value p " pp 1 , p 2 q :" EtP u of the random variable P is identified, in a classical manner [START_REF] Ogden | Fitting hyperelastic models to experimental data[END_REF], by fitting the mean experimental response with an Ogden-type constitutive model of orders m " n " 1: pp, ηq :" argmin pa,bq P AˆB J pa, bq ,

where the cost function J is given by J pa, bq :"

N ÿ k"1 ´σpυ exp k ; a, bq ´ξexp k ¯2 , (37)
and σpυ exp k ; a, bq is the predicted mean stress at the stretch υ exp k , obtained by substituting pa, bq for pp, ηq in Eq. ( 8) (with m " n " 1). The solution search is performed over the admissibles sets

A " pR ˚q2 (38)
and

B " tη P R 2 , η 1 ě 2 , η 2 ě 3{2u (39)
induced by Eqs. (9a)-(9c). This optimization problem can be solved by using a standard algorithm for constrained nonlinear problems.

Defining

Step 2: calibration of pρ 1 , ρ 2 q and pλ 1 , λ 2 q 270 Let us address the calibration of the remaining hyperparameters pρ 1 , ρ 2 q and pλ 1 , λ 2 q. From Eq. ( 28) and the statistical independence of µ and P j , we deduce that

p 1 " 2Etµu η 2 1 EtU u ( 40
)
and

p 2 " 2Etµu η 2 2 p1 ´EtU uq . ( 41
)
Since µ " Gpρ 1 , ρ 2 q and U " Dirpλ 1 , λ 2 q (note that the Dirichlet distribution corresponds to the Beta distribution for univariate random variables), one has

µ " ρ 1 ρ 2 (42)
and

p 1 " 2µ η 2 1 λ 1 λ 1 `λ2 , p 2 " 2µ η 2 2 λ 2 λ 1 `λ2 . ( 43
)
Solving for ρ 2 and λ 1 in Eq. ( 42) and in the first equation of Eq. ( 43), it is found that

ρ 2 " µ ρ 1 , (44)
and

λ 1 " η 2 1 p 1 η 2 2 p 2 λ 2 " u 1 ´u λ 2 , ( 45
)
where u :" EtU u and use was made of Eqs. ( 40)- [START_REF]Stochastic modeling of the ogden class of stored energy functions for hyperelastic materials: the compressible case[END_REF]. The idea in this final step is to define the parameters ρ 2 and λ 1 through Eqs. ( 44) and (45) (recall that tη k u 2 k"1 and tp k u 2 k"1 are known from Step 1), and to consider ρ 1 and λ 2 as the only hyperparameters to be calibrated. For this purpose, let us introduce the cost function K defined as Kpx, yq :"

N ÿ k"1 pς Σ pυ exp k ; x, yq ´ςexp k q 2 , ( 46
)
where ς Σ pυ exp k ; x, yq is the standard deviation of the stress Σ at stretch υ exp k (which can be estimated by using the closed-form expressions provided at the end of section 3.2.1), with

µ " G ´x, µ x ¯, (47)
and

U " Dir ˜η2 1 p 1 η 2 2 p 2 y, y ¸. (48)
The hyperparameters ρ 1 and λ 2 are then defined as the solution of the minimization problem: pρ 1 , λ 2 q :" argmin px,yq P s0,`8rˆs1,`8r

Kpx, yq .

In this work, the optimization problems defined by Eq. ( 36) and Eq. ( 49) are solved by using the Matlab function fmincon (with an sqp algorithm) for a set of 1, 000 initial guesses determined by Latin Hypercube sampling. Here, the randomization of the starting points aims at reducing 275 in part the impact of the local optimization technique.

Discussion

A number of remarks regarding the calibration strategy are relevant at this point. First of all, the generalization of the calibration method to arbitrary orders shows that 280 the dimension of the problem related to the calibration of ρ 1 and λ np (which corresponds to λ 2 in the previous section) is, indeed, independent of the values of m and n. As a consequence, the methodology does not suffer from a curse of dimensionality and turns out to be very robust in 285 that sense.

Secondly, and while the set of model parameters can be shown to be unique when the constraint equations (given by Eqs. ( 14), ( 15), ( 20) and ( 21)) can be all estimated from experiments, considering the mean and the variability of the stress response does not allow the hyperparameters to be uniquely defined. In fact, it should be observed that the random stress involves nonlinear functions of µ and U . Therefore, a given variability in the random stress may equivalently be generated by different combinations of fluctuations in µ and U (see Eq. ( 32), ( 33), ( 34) and ( 35)). Note that issues related to (non)uniqueness are also encountered in a deterministic context (which is relevant to the identification of pη 1 , η 2 q and pp 1 , p 2 q within Step 1); see the discussion in [START_REF] Ogden | Fitting hyperelastic models to experimental data[END_REF] for instance.

Additionally, the choice of the hyperparameters to be identified within Step 2 turns out to be arbitrary, so that the associated optimization problem may alternatively be formulated by selecting any couple pρ k , λ q for 1 ď k ď 2 and 1 ď ď 2.

Finally, it should be noted that Eq. ( 45) can be used in order to investigate the sensitivity of the cost function K (see Eq. ( 46)) with respect to the second design variable. Indeed, the coefficient of variation δ U of U is given by

δ U " d λ 2 λ 1 pλ 1 `λ2 `1q , (50)
with

λ 1 " u 1 ´u λ 2 . (51)
When the mean value u (which can be estimated within

Step 1 by using Eq. ( 40) for example) takes small positive values (it is recalled that 0 ă U ă 1 almost surely), one has λ 1 ! λ 2 and U exhibits large statistical fluctuations. Similarly, it is seen that λ 1 " λ 2 when u is close to 1, in which case the fluctuations of U may become very small. In both cases, the coefficient of variation δ U quickly reaches a plateau as λ 2 becomes larger, and this behavior may decrease the robustness of the formulation within Step 2.

Applications to soft biological tissues

In this section, the stochastic stored energy functions (with m " n " 1) and the identification methodology are applied in order to model the variability exhibited by various soft biological tissues, including brain and liver tissues, as well as spinal cord white matter.

Brain tissue

In this first application, we address the modeling of the variability exhibited by brain tissues. Toward this aim, we consider the experimental results provided in [START_REF] Jin | A comprehensive experimental study on material properties of human brain tissue[END_REF]. The database is composed of 72 samples of brain tissues tested for strains ranging up to 50% and for strain rates 9 ε P t0.5, 5, 30u s ´1. Below, we address the calibration task in unconfined compression for low and medium strain rates. Experimental data are given for white and gray matter, and final data are expressed in terms of mean and standard deviation plots at increasing stretches.

The optimal mean parameters and exponents pp, ηq, computed within Step 1 (see section 3.2), are reported in Tables 1 and2

335

It is seen that the inequality constraints raised by the polyconvexity and coerciveness properties (see Eqs. (9a)-(9c)) are fulfilled, so that the mean model is admissible. Regarding the computation of the Lagrange multipliers ρ 1 and λ 2 , performed within the second step of the methodol-340 ogy, the algorithm is found to converge in a few iterations, regardless of the initial guess. The mean and variability of the random stress predicted by the calibrated stochastic model are quantitatively compared with experimental data in Figs. 1, 2, 3 and4, where the confidence interval 345 at 95% is also reported.

It can be observed that the first-order stochastic stored energy function can reproduce the experimental results very well for both gray and white matters at strain rates 9 ε P t0.5, 5, 30u s ´1. Black lines: experimental data extracted from [START_REF] Jin | A comprehensive experimental study on material properties of human brain tissue[END_REF]. Red lines: calibrated stochastic model. Yellow area: 95% confidence region.

Liver tissue

As a second example, we now turn to the modeling of bovine liver tissue. Specifically, the statistical results are extracted from [START_REF] Roan | The nonlinear material proper-515 ties of liver tissue determined from no-slip uniaxial compression experiments[END_REF], in which uniaxial compression tests were performed at a 0.01 s ´1 strain rate. It should be noted that the (deterministic) fitting in the above reference was obtained by using an exponential-type stored energy function. The calibrated parameters are listed in Table 3. A qualitative comparison between the statistical properties estimated by means of Monte Carlo simulations with the calibrated stochastic model and those computed from the experiments can be visualized in Fig. 5. Similarly to the case of brain tissues, it is seen that the probabilistic model and the identification strategy allows the data to be reproduced with a reasonably high level of accuracy, hence 365 showing the relevance of the overall methodology.

Spinal cord white matter

In this final application, we address the modeling of the variability exhibited by porcine spinal cord white matter under unconfined compression. Experimental data can be 370 found in [START_REF] Sparrey | Compression behavior of porcine spinal cord white matter[END_REF], where the mechanics of spinal cord injury and its finite element implementation are investigated. The database is composed of 104 independent samples extracted from Yorkshire pigs. Unconfined compression tests were performed for strains ranging up to 40% and for a strain

Experimental data Stochastic model Confidence region

Figure 5: Graph of υ Þ Ñ Σpυq (kPa) and υ Þ Ñ ς Σ pυq (kPa) for liver tissues in uniaxial compression at 9 ε " 0.01 s ´1. Black lines: experimental data extracted from [START_REF]Strain rate-dependent viscohyperelastic constitutive modeling of bovine liver tissue[END_REF]. Red lines: calibrated stochastic model. Yellow area: 95% confidence region.

are listed in Table 4 and the predicted results are compared with the experimental ones in Figs. 6, 7 and 8.

A very good agreement is again observed between the predictions of the stochastic model and the experimental results.

Conclusion

In this paper, we have addressed the construction and inverse identification of stochastic constitutive models for hyperelastic biological tissues. The approach builds upon information theory, which essentially allows the constraints related to existence theorems in nonlinear elasticity to be taken into account. This ingredient is key to deriving a stochastic constitutive law that is both physically sound and mathematically consistent. The model can, indeed, be seen as a stochastic extension of the Ogden class of polyconvex stored energy functions, and relies on a low-

  105

Figure 1 :

 1 Figure 1: Graph of υ Þ Ñ Σpυq (kPa) and υ Þ Ñ ς Σ pυq (kPa) for gray matter tissues in unconfined compression at 9ε " 0.5 s ´1. Black lines: experimental data extracted from[START_REF] Jin | A comprehensive experimental study on material properties of human brain tissue[END_REF]. Red lines: calibrated stochastic model. Yellow area: 95% confidence region.

Figure 2 :

 2 Figure2: Graph of υ Þ Ñ Σpυq (kPa) and υ Þ Ñ ς Σ pυq (kPa) for gray matter tissues in unconfined compression at 9 ε " 5 s ´1. Black lines: experimental data extracted from[START_REF] Jin | A comprehensive experimental study on material properties of human brain tissue[END_REF]. Red lines: calibrated stochastic model. Yellow area: 95% confidence region.

Figure 3 :

 3 Figure 3: Graph of υ Þ Ñ Σpυq (kPa) and υ Þ Ñ ς Σ pυq (kPa) for white matter tissues in unconfined compression at 9ε " 0.5 s ´1. Black lines: experimental data extracted from[START_REF] Jin | A comprehensive experimental study on material properties of human brain tissue[END_REF]. Red lines: calibrated stochastic model. Yellow area: 95% confidence region.

Figure 4 :

 4 Figure4: Graph of υ Þ Ñ Σpυq (kPa) and υ Þ Ñ ς Σ pυq (kPa) for white matter tissues in unconfined compression at 9 ε " 5 s ´1. Black lines: experimental data extracted from[START_REF] Jin | A comprehensive experimental study on material properties of human brain tissue[END_REF]. Red lines: calibrated stochastic model. Yellow area: 95% confidence region.

Figure 6 :

 6 Figure 6: Graph of υ Þ Ñ Σpυq (kPa) and υ Þ Ñ ς Σ pυq (kPa) for spinal cord white matter in unconfined compression at 9ε " 0.005 s ´1. Black lines: experimental data extracted from[START_REF] Sparrey | Compression behavior of porcine spinal cord white matter[END_REF]. Red lines: calibrated stochastic model. Yellow area: 95% confidence region.

Figure 7 :

 7 Figure 7: Graph of υ Þ Ñ Σpυq (kPa) and υ Þ Ñ ς Σ pυq (kPa) for spinal cord white matter in unconfined compression at 9ε " 0.05 s ´1. Black lines: experimental data extracted from[START_REF] Sparrey | Compression behavior of porcine spinal cord white matter[END_REF]. Red lines: calibrated stochastic model. Yellow area: 95% confidence region.

Figure 8 :

 8 Figure 8: Graph of υ Þ Ñ Σpυq (kPa) and υ Þ Ñ ς Σ pυq (kPa) for spinal cord white matter in unconfined compression at 9ε " 0.5 s ´1. Black lines: experimental data extracted from[START_REF] Sparrey | Compression behavior of porcine spinal cord white matter[END_REF]. Red lines: calibrated stochastic model. Yellow area: 95% confidence region.

  A, A and rAs. The table of main notations is reported below (here, stresses are defined for a uniaxial compression test).

			Σ	Random principal first Piola-Kirchoff stress
			Et¨u	Mathematical expectation
			}rAs}	Frobenius norm of rAs
		90	Finally, k 0 denotes a normalization constant, the value of
	w	Deterministic stored energy function
	W	Random stored energy function
	µ	Deterministic shear modulus
	µ	Random shear modulus
	σ	Principal first Piola-Kirchoff stress

Table 1 :

 1 for the compression tests at different strain rates. calibrated parameters (gray matter in compression).

	9 ε	0.5 s ´1	5 s ´1
	p 1 (kPa)	0.1467	0.1881
	p 2 (kPa)	0.0457	0.7823
	η 1	5.5945	4.5907
	η 2	1.991	1.5
	µ (kPa)	2.3863	2.3863
	λ 1	253.5375 22.461
	λ 2	9.9982	9.9732
	ρ 1	2.3679	6.6398
	ρ 2	1.0078	0.4311
	9 ε	0.5 s ´1	5 s ´1
	p 1 (kPa) 1.1131	1.4720
	p 2 (kPa) 1.1120	0.4897
	η 1	2	2.6185
	η 2	1.5	1.5
	µ (kPa)	3.4772	5.5973
	λ 1	3.2290 91.5357
	λ 2	1.8146	9.9919
	ρ 1	33.2721 6.6819
	ρ 2	0.1045	0.8377

Table 2 :

 2 calibrated parameters (white matter in compression).

Table 3 :

 3 calibrated parameters (liver tissues in compression).

	375
	rate 9 ε P t0.005, 0.05, 0.5u s ´1. The calibrated parameters

Table 4 :

 4 calibrated parameters (spinal white matter in compression).

	9 ε	0.005 s ´1	0.05 s ´1 0.5 s ´1
	p 1 (kPa) 1.0003 ˆ10 ´10	0.0102	0.0019
	p 2 (kPa)	0.0518	0.0060	0.0111
	η 1	2.1306	10.3693 11.4817
	η 2	4.5327	4.2260	5.3651
	µ (kPa)	0.5317	0.6048	0.2862
	λ 1	4.2683 ˆ10 ´9	10.2379	7.9405
	λ 2	9.9991	1	10
	ρ 1	25.3357	16.6017 23.7158
	ρ 2	0.0210	0.0364	0.0121

Appendix A.1. General derivations

Let A be a vector-valued random variable defined by a probability density function p A . We assume that the probability density function p A , defined from R n into R `, has a support which is denoted by S A (hence, p A paq " 0 for a R S A ), with S A Ď R n . It is assumed that some objective information is available and can be formalized in the form of a mathematical expectation:

where h is a given measurable mapping from R n into R N , with N ě 1, and f is a given vector in R N . The MaxEnt principle then states that the most objective model, given the above constraints, is given by [START_REF]Information theory and statistical mechanics i[END_REF][START_REF] Jaynes | Information theory and stastitical mechanics ii[END_REF][START_REF] Shannon | A mathematical theory of communication[END_REF] p A " argmax

where Stpu :" ´żR n ppAq log pppAqq da (A. [START_REF] Balzani | A polyconvex framework for soft biological tissues. adjustment 425 to experimental data[END_REF] is the so-called Shannon's entropy of p and C A is the set of all probability density functions, with support S A , satisfying the constraints defined by Eq. (A.1). In order to solve the above functional optimization problem, a Lagrange multiplier τ with values in an admissible set T Ă R N and associated with Eq. (A.1) is introduced. Proceeding with the calculus of variation on the associated Lagrangian (see e.g. [START_REF] Luenberger | Optimization by vector space methods[END_REF]), the solution is found to be defined as

where 1 S A is the indicator function of S A and k 0 is a normalization constant such that ş S A p A paq da " 1.

405

Appendix A.2. Application 1: probabilistic modeling of the random shear modulus µ

It is assumed that the constraints related to the random shear modulus µ are given by Eqs. ( 15) and ( 16). Hence, one has

with n " 1 and N " 2 in this particular case. Let p1 τ1 q and τ 2 be two Lagrange multipliers associated with the two constraints. Using Eq. (A.4), it is found that the probability density function p µ writes p µ pµq " 1 R ˚pµq k 0 exp p´τ 2 µ ´p1 ´τ1 q logpµqq . (A.6) Upon using the changes of variables ρ 1 " τ 1 and ρ 2 " 1{τ 2 , it is seen that p µ can be written as

As a result, the random shear modulus µ is a Gamma random variable with parameters ρ 1 ą 0 and ρ 2 ą 0.

Appendix A.3. Application 2: probabilistic modeling of the 410 normalized random variable U

The probability density function p U of random variable U is constructed by considering that the only available information is given by Eqs. ( 20) and ( 21), with n " n p ´1, N " n p and

Let p1 ´λk q, for 1 ď k ď n p , be n p Lagrange multipliers associated with the above constraints. Applying the Max-Ent principle (see Eq. (A.4)), it is found that p U takes the form

where the support S is given by Eq. ( 19). After little algebra, it is found that p U reads as

(A.10) As a result, see e.g. [START_REF] Johnson | Continuous 470 multivariate distributions, volume 1, models and applications[END_REF], random variable U follows a Dirichlet distribution with parameters pλ 1 , λ 2 , . . . , λ np q. Note that the Dirichlet distribution corresponds to a classical Beta distribution for n p " 2.


Francesco Viola

P.-T Brun

Benjamin Dollet

François Gallaire

Foam on troubled water: capillary induced finite-time arrest of sloshing waves

Interfacial forces exceed gravitational forces on a scale small relative to the capillary length -two millimeters in the case of an air-water interface -and therefore dominate the physics of sub-millimetric systems. They are of paramount importance for various biological taxa and engineering processes where the motion of a liquid meniscus induces a viscous frictional force that exhibits a sublinear dependence in the meniscus velocity, i.e. a power law with an exponent smaller than one. Interested in the fundamental implications of this dependence, we use a liquid-foam sloshing system as a prototype to exacerbate the effect of sublinear friction on the macroscopic mechanics of multi-phase flows. In contrast to classical theory, we uncover the existence of a finite-time singularity in our system yielding the arrest of the fluid's oscillations. We propose a minimal theoretical framework to capture this effect, thereby amending the paradigmatic damped harmonic oscillator model. Our results suggest that, although often not considered at the macroscale, sublinear capillary forces govern the friction at liquid-solid and liquid-liquid interfaces.

I. INTRODUCTION

No one who has glimpsed at Hokusai's iconic wave paintings can fail to notice the truthful yet stylized rendering of the singular nature of such multi-scale flows. On a theoretical point of view, overturning, breaking and atomization are supported by nonlinear terms in the equations prescribing the waves' dynamics at large amplitude, so that finding their exact solution is often arduous. Nonlinearities, however, usually become of negligible relative importance for sufficiently low amplitudes of motion, thereby enabling perturbative linearized approaches [START_REF] Lamb | Hydrodynamics[END_REF] . These theoretical linear problems are easier to solve than their nonlinear counterparts and have been extremely successful at predicting the stability, characteristic frequencies and growth rates of a wide range of physical systems arising in natural and industrial settings, a classic example of which is the sloshing of fluid in a container [START_REF] Ibrahim | Liquid Sloshing Dynamics: Theory and Applications[END_REF] . In such a case, the excess of gravitational potential energy associated to the deformation of the free surface is progressively relaxed when periodically turned into kinetic energy and damped by viscous forces [START_REF] Ibrahim | Liquid Sloshing Dynamics: Theory and Applications[END_REF][START_REF] Case | Damping of surface waves in an incompressible liquid[END_REF] . In the small amplitude regime, the fundamental sloshing frequency for a cylindrical container is known to be ω 0 = 1.84 g/R where g is the acceleration of gravity and R the container radius [START_REF] Ibrahim | Liquid Sloshing Dynamics: Theory and Applications[END_REF] . The classical contributions to the viscous damping rate, σ, have been identified as the dissipation in the oscillating boundary layers near the container walls and the free surface, and the bulk dissipation [START_REF] Case | Damping of surface waves in an incompressible liquid[END_REF] . Those dissipative forces have different scalings with respect to the governing parameters (radius R, viscosity µ, density ρ, gravity g), but all depend linearly on the free surface velocity [START_REF] Ibrahim | Liquid Sloshing Dynamics: Theory and Applications[END_REF][START_REF] Case | Damping of surface waves in an incompressible liquid[END_REF] .

Here, we study the effect of a small modification to this classical problem, as we place a thin layer of foam on a volume of water (see Fig. 1), and explore the dramatic change a) Electronic mail: francois.gallaire@epfl.ch in its mechanistic properties. It has been recently shown [START_REF] Sauret | Damping of liquid sloshing by foams[END_REF] that sloshing in a rectangular basin could be strongly damped by the addition of a shallow foam layer at the top. While foam is a complex media with visco-elasto-plastic rheological properties 5 , this significant attenuation was attributed to friction between the foam and the wall, which dominates the dissipation in the foam bulk [START_REF] Sauret | Damping of liquid sloshing by foams[END_REF] . In this work, we focus our attention on the nonlinear nature of the dissipation occurring in the foam layer. This source of dissipation emanates from the relative displacement of the liquid channels separating the foam bubbles, known as Plateau borders [START_REF] Cantat | Liquid meniscus friction on a wet plate: Bubbles, lamellae, and foams[END_REF] . The motion of these multiple air-liquid interfaces along the container results in a frictional force that scales nonlinearly with their sliding velocity. Indeed, this force scales like Ca 2/3 , where Ca is the capillary number based on the sliding velocity [START_REF] Cantat | Liquid meniscus friction on a wet plate: Bubbles, lamellae, and foams[END_REF][START_REF] Bretherton | The motion of long bubbles in tubes[END_REF][START_REF] Denkov | Wall slip and viscous dissipation in sheared foams: Effect of surface mobility[END_REF] . Owing to their sublinearity (the exponent on the velocity is smaller than one), these frictional terms are found to dominate all linear contributions for small enough values of Ca, which naturally arise for small amplitudes of motion. Therefore, unlike nonlinearities with exponents greater than one, they cannot be neglected, thereby ruining all linear theoretical approaches. Unlike the picture of a linearly damped oscillator classically used in textbooks to model sloshing, we show experimentally and theoretically that the sloshing of the liquid phase does not relax exponentially. The damping rate in fact increases catastrophically at small amplitude and yields a finite-time arrest of the fluid's motion.

II. EXPERIMENTS

The setup consists of a cylindrical Plexiglas container of internal radius R = 37.05 mm that is sealed to an aluminum plate connected itself to a hinge. The container is filled with a column of water of height H which is surmounted by a foam layer of thickness h, see Figure 1A. The foam is produced when injecting a commercial hand soap (Hypron SA) through a nozzle with a grid with pore size 0.2 mm and has typical bubbles diameter d = 0.25 mm. The foam's liquid fraction is measured to be 25% once poured on top of the water column (evaluated using the weight of a volume of foam). This large value relative to dry foams, is due to capillary forces yielding an intake of water effectively enlarging the foam Plateau borders 5 . In the following, the liquid fraction across the foam layer, whose thickness ranges from h = 2.5 mm to 6.5 mm, is assumed uniform. Using an oscillating drop tensiometer (Tracker, Teclis), we measured the surface tension, γ = 31.5 mN/m, of the soap solution. Its viscoelastic surface dilatational modulus, quantifying the relative variation of surface tension upon dilatation and compression, is found below resolution (0.1 mN/m), showing that the interface can be considered as mobile [START_REF] Denkov | Wall slip and viscous dissipation in sheared foams: Effect of surface mobility[END_REF] and that free shear boundary conditions apply.

Experiments are initiated by tilting the container of an angle ψ and then releasing it, see Fig. 1A. The dynamics following the impact is recorded with a high-speed camera (Phantom Micro M310) operating at 300 Hz. In order to investigate the waves dynamics at small scales, we coupled a Nikon 50 mm lens with an optical bellows obtaining a resolution of 2•10 per pixel on the water-foam interface displacement. After a short transient regime, the flow response was found to exhibit a single mode of oscillation (see Fig. 1B). The frequency of oscillations, 3.50 Hz, matches the theoretical inviscid prediction ω 0 = 1.84g/R within 1%. This result was expected for h = 0 as the deep water condition H > 2R is satisfied and the container size is well above the capillary length 1 , i.e. γ/ρgR [START_REF] Ibrahim | Liquid Sloshing Dynamics: Theory and Applications[END_REF] 1. This frequency is found unchanged when adding foam (0 < h < 6.5 mm), as already observed in previous experiments [START_REF] Sauret | Damping of liquid sloshing by foams[END_REF] and discussed in section III B.

The attenuation of the gravity waves is investigated by recording the motion X(t) of the water-foam interface at the azimuth corresponding to the largest wave oscillation. Using the function 'orthogonal view' in ImageJ, we follow the evolution of a pixel line over time, thereby reconstituting the oscillatory motion X(t) of the contact line at a given location, see Fig. 1B. The natural logarithm D = ln(X n /X 0 ) of the ratio between the amplitude of the nth peak, X n , and the amplitude of the first peak, systematically chosen to X 0 = 5 mm, is used as an instantaneous measure of the growth rate. For the case h = 0 (no foam), we recover the classical result of a constant damping rate. The logarithmic decrement is linear in time for most of the dynamics, D ∝ -σω 0 t, as indicated in Fig. 2A (see section IV for a rationalization of the shaded area). We found σ = 0.0098, which exceeds the theoretical expression [START_REF] Case | Damping of surface waves in an incompressible liquid[END_REF] of the viscous damping rate by a factor 2. This result is in line with other experimental measurements in the literature [START_REF] Cocciaro | Capillarity effects on surface gravity waves in a cylindrical container: wetting boundary conditions[END_REF] and is commonly attributed to the advancing and receding of the dynamic contact line 2 . As soon as a thin foam layer is added in the system, D is no longer linear, as shown in Fig. 2A for values of h ranging from 2.5 mm to 6.5 mm. We find D to be concave down and negatively correlated to h. In other words, the damping rate is greater with foam as evidenced by previous experiments 4 but also evidently increases as the amplitude of oscillations decreases. In particular this result translates into the fact that the interface does not relax exponentially at a given rate, but instead its decay is progressively hastened. This novel observation suggests the existence of a mechanism submitting smaller oscillations to a stronger damping, that we aim to rationalize.

III. A NONLINEAR FLUID OSCILLATOR A. Equivalent mechanical system

From this point onward, we restrict ourselves to thin layers of foam and we assume the entire layer of height h is oscillating along with the liquid phase. We model the total oscillating mass as πR 2 (ρλ+ρ f h), where λ ∼ R is the penetration length of the wave 1 and ρ and ρ f denote the density of water and foam, respectively. Given the physical parameters of the problem we place ourselves in the limit ρ f h ρλ, i.e. we neglect foam's inertia relative to the one of water. Similarly, the restoring gravitational force acting on this volume of fluid may be estimated considering solely the liquid phase: F g ∼ ρgπR 2 X, where X is the amplitude of the water-foam interface and g the acceleration of gravity. In We now turn to model the force F f associated with capillary effects in the foam layer. Let us consider a single wall Plateau border, i.e. the edge separating neighboring bubbles at the wall, sliding along the container vertical wall with velocity Ẋ. Damping of its motion originates from the transition region between the wetting film at the wall and the moving wall Plateau border [START_REF] Cantat | Liquid meniscus friction on a wet plate: Bubbles, lamellae, and foams[END_REF] . The matching between the outer meniscus and the transition zone, involving a balance between viscous and capillary forces, defines the friction force per unit length of a wall Plateau border scaling as γCa 2/3 , where Ca = µ Ẋ/γ is the capillary number and γ denotes the surface tension [START_REF] Cantat | Liquid meniscus friction on a wet plate: Bubbles, lamellae, and foams[END_REF] . The total length of surface Plateau borders is given by their typical length, d, times the number of Plateau borders in contact with the container that scales as 2πRh/d 2 . Thus, the friction force acting at the wall in the foam layer reads

F f ∼ h d Rγ 1/3 µ 2/3 Ẋ| Ẋ| -1/3 , (1)
where the absolute value is introduced to conveniently treat negative speeds. This result is reminiscent to the force acting on a plate pulled out of a bath [START_REF] Bretherton | The motion of long bubbles in tubes[END_REF][START_REF] Landau | Dragging of a liquid by a moving plate[END_REF] and is valid for a mobile interface, i.e. an interface where a free shear boundary condition applies. This point has been verified for the foam used in our experiments (see section II).

Using the estimation for the volume of fluid effectively oscillating and the forces acting on it, we write Newton's 2 nd law:

πR 2 (ρλ) Ẍ = F g + F w + F f . (2)
Equation ( 2) constitutes a nonlinear oscillator where inertia is resisted by the restoring force of gravity and the two dissipative terms, F w and F f (see the equivalent mechanical system depicted in Figure 3). We note that F w varies linearly with Ẋ, while F f has a sublinear dependence on the interface velocity Ẋ. As a consequence, F f is expected to dominate the dynamics when Ẋ is sufficiently small. After a renormalization by the inertial term, eq. ( 2)

yields Ẍ + ω 2 0 X = -2σω 0 Ẋ -α ω 4/3 0 R 1/3 Ẋ| Ẋ| -1/3 , (3)
where ω 0 is the angular frequency scaling as ω 0 ∼ g/R, and σ is the viscous damping coefficient which is proportional to σ ∼ √ ν/(R 3/4 g 1/4 ). These scaling arguments are in line with the theoretical results obtained using first principles [START_REF] Lamb | Hydrodynamics[END_REF][START_REF] Case | Damping of surface waves in an incompressible liquid[END_REF] . The last term on the right hand side is the nonlinear damping associated with the foam layer, where the prefactor ω 4/3 0 R 1/3 results from the renormalization, and α is the dimensionless coefficient of sublinear damping, which scales as α ∼ hd -1 γ 1/3 µ 2/3 ρ -1 R -5/3 g -2/3 .

B. Asymptotic solution and comparison to experiments

In the limit of negligible capillary effects, α = 0, equation ( 3) reduces to a linearly damped harmonic oscillator. This equation is traditionally used to model the free-surface oscillations of a liquid in a tank 2 , with solution:

X(t) = X 0 e -ω0σt cos( 1 -σ 2 ω 0 t), (4)
where X 0 is the initial amplitude. This expression is consistent with the experimental data reported in Fig. 2 for h = 0 (no foam). With foam and thus α = 0, eq. ( 3) does not have a closed-form analytical solution. However, an asymptotic solution can be obtained by using the method of multiple scales as detailed in Nayfeh [START_REF] Nayfeh | Perturbation Methods[END_REF] . In the limit of small damping, σ, α 1, the solution of eq. ( 3) is asymptotic to X(t) ∼ a(t) cos(ω 0 t), where the sinusoidal term is the solution of the undamped oscillator (σ = α = 0 as in the leading order problem in the asymptotic scheme). The slow varying envelope a(t) is determined at the next order by imposing a compatibility condition that involves only the right hand side terms in eq. ( 3) oscillating at the natural frequency of the system ω 0 . X(t) reads:

X(t) = X 1/3 0 - α σ cR 1/3 e ω0σt/3 -1 3 e -σω0t a(t)
cos(ω 0 t).

(

) 5
See supplemental material 11 at [URL] for a detailed derivation of the asymptotic solution and the definition of the constant c that is the Fourier coefficient of the Fourier component of Ẋ| Ẋ| -1/3 synchronized to the natural frequency of the system ω 0 . Note that the frequency shift √ 1 -σ 2 in eq. ( 4), is not present in the asymptotic solution (5) as it belongs to a higher order.

The asymptotic solution (5) to the sublinear model (SLM), is now compared and fitted to our experimental results. In Figure 2(b) the coefficients ω 0 and σ are set to their reference values obtained experimentally in the h = 0 limit (ω 0 = 22.00 rad/s, σ = 9.8 • 10 -3 ). The sublinear coefficient α is thus the only unknown parameter of the model. Here, α is determined independently for each of the h = 0 cases by fitting the logarithmic decrement with equation (5). We find that our model captures the concave aspect of the linear decrement measured from experiments. Additionally, the larger the film thickness, the stronger the nonlinear damping coefficient. More precisely, our data in Fig. 2B shows that the damping coefficient α increases linearly with the foam thickness h, further validating the proposed physical picture: in the range of thin foam layers investigated, the foam acts as collection of Plateau borders, which all equally contribute to the effective dissipative force. The number of contact lines, earlier estimated as the foam thickness to bubble diameter ratio h/d, is thus the relevant scale to evaluate α.

C. Sublinear friction and its finite-time arrest

Assured that the reduced model (3) and its solution (5) are a quantitative representation of the sloshing dynamics with foam, we turn to examine the mechanistic impact of 6) with α(h) reported in Fig. 2B). Symbols correspond to the values deduced from the experiments.

the sublinear damping introduced in our model. In Fig. 4A we report the sublinear model (SLM) envelope a(t) in the particular, yet representative, case h = 4 mm. Its concave down behavior contrasts with the classical exponential relaxation adopted by previous works [START_REF] Sauret | Damping of liquid sloshing by foams[END_REF] , indicated by the blue line. This effect is due to an increase of the relative importance of sublinear friction forces compared to the inertial and restoring terms as the amplitude of oscillation decreases. The dominance of interfacial effects is induced by their sublinear dependance in the oscillating velocity. In fact, they diverge in the limit of vanishing oscillations, yielding a finite-time singularity. For the sake of the validation of our asymptotic solution (5), we report also in Fig. 4A the numerical solution of eq. ( 3) obtained using a third order Runge-Kutta scheme. Upon examination of the behavior of eq. ( 5) at very small amplitudes, we find the amplitude of sloshing to be exactly nil for time t = t * , where t * is the time of arrest. This behavior contrasts the classical linear theory giving rise to an exponential decay for which t * = ∞. In our model, t * is finite and using eq. ( 5), we can derive its analytical expression:

t * = 3 σω 0 ln 1 + c -1 σ α(h) X 0 R 1/3 . (6)
In Fig. 4B we show the variation of t * with h and report the values corresponding to our experiments (as deduced from the fit of our experimental data using the SLM). The time of arrest t * diverges in the limit where h = 0, as expected since this case corresponds to an exponential relaxation in our theoretical framework (see eq. ( 4)). A small foam layer (h = 2.5 mm) is sufficient to bring this time to a finite and relatively small value (t * 12 s). Other values of t * corresponding to larger values of h are of the same order (the smallest is 5 s) and decrease if the foam thickness h increases.

IV. CONCLUSIONS

In our experiments we found that adding foam to a liquid dramatically damps its oscillatory motion but hardly affects its sloshing frequency. Without foam, sloshing is well approximated by a linearly damped harmonic oscillator relaxing exponentially with characteristic time t µ = 1/(σω 0 ), thereby predicting an arrest of motion for t * = ∞. In contrast, owing to its cellular structure, foam introduces sublinear capillary forces, which, in turn, lead to a finite-time arrest of the fluid's motion. Equation ( 6) provides the expression for the time of arrest t * , highlighting the importance of the ratio between the viscous damping σ ∼ ν 1/2 g -1/4 R -3/4 and the sublinear damping α ∼ hd -1 γ 1/3 µ 2/3 R -5/3 g -2/3 ρ -1 coefficients. The time of arrest t * is negatively correlated to α/σ and diverges when this ratio is small (corresponding to a linearly damped oscillator). Specifically, an increase in h, an increase in γ or a decrease in d, all yield a decrease in t * . The scaling argument for the sublinear coefficient α is in good agreement with the experimental measurements for various foam height h, see section III B. On the other hand, its dependence on the fluid density, ρ, its viscosity, µ, and surface tension, γ, along with the dependence on the container size, R, and the foam bubble diameter, d, can in principle be probed experimentally by using different fluids, container size and foams and it is a future perspective of this work.

Note that this expression is only valid for small but finite foam layers. Sloshing with a large amount of foam may lead to the situation where only a fraction of the total volume of foam is effectively in motion. The upper part of the foam layer remains still, hence does not participate to the sublinear damping examined in this work. We therefore anticipate that α(h) does not remain linear for arbitrarily large values of h, but instead saturates, in turn giving rise to a plateau for t * in the limit of large foam thicknesses.

Similarly, the case d = 0 and h = 0 is ill-defined in our model. For these values, corresponding to a simple liquid-air interface, a liquid meniscus forms at the container walls. During sloshing, this advancing and receding contact line will generate a small, yet nonlinear, friction force 5 , which affects the damping rate measured in the experiments [START_REF] Ibrahim | Liquid Sloshing Dynamics: Theory and Applications[END_REF][START_REF] Cocciaro | Capillarity effects on surface gravity waves in a cylindrical container: wetting boundary conditions[END_REF] . Evidence of these effects is apparent at very small amplitude as highlighted by the shaded area in Fig. 2A for pure water. Although the accuracy of the measurements does not allow for a systematic analysis at this stage, this observation suggests that our work could generalize to a single moving contact line. Accounting for those effects -without foam to exacerbate them -requires a realistic contact line model 2 , with features such as contact angle hysteresis and the variations of contact angle with the contact line velocity.

Finally, we note that finite-time arrests are commonly observed in mechanical systems involving solid friction as a consequence of Coulomb's law, but have not frequently been reported in fluidic systems. Their importance in practical situations such as Euler's disk has been the subject of investigations in recent years [START_REF] Moffatt | Euler's disk and its finite-time singularity[END_REF] . Here, we have demonstrated that the effect of pure fluid dissipation in sub-millimetric menisci at a contact line is sufficient to bring a body of fluid to rest in finite time albeit the fluid's volume is significantly larger than that of the meniscus.

  FIG. 1. (A) Experimental setup: a cylindrical container is tilted by an angle ψ to initiate the sloshing of a volume of water surmounted by foam. The foam layer of thickness h H and h R comprises bubbles of typical radius d. (B) Shown is the spatio-temporal view of the water-foam contact line X(t) captured in a typical experiment.

1 FIG. 2 .

 12 FIG. 2. (A)Logarithmic decrement, D = ln(Xn/X0), of the sloshing amplitude for foam thicknesses ranging across values h = 0, 2.5, 3.3, 4, 5, 6.5 mm. Symbols indicate measurements, while continuous lines represent the envelope solutions of the sublinear model (SLM), with a single fitting parameter α (see a(t) in eq. 5) . (B) Shown is the dependance of the friction force exerted by the foam layer, F f , through the friction coefficient α as found in experiments. F f increases linearly with the foam thickness h, i.e. with the number of layers of oscillating bubbles.

FIG. 3 .

 3 FIG. 3. A sketch of our reduced model for sloshing. The effects of foam are accounted by the sublinear term F f ∝ Ẋ| Ẋ| -1/3 . This term amends the classical damped oscillator model where the inertial term is balanced by the restoring force, Fg ∝ X, and the viscous damping, Fw ∝ Ẋ.

  absence of foam, according to the seminal work of Case & Parkinson 3 , three main dissipation sources exist, they acting: (i) at the free surface, (ii) in the fluid bulk and (iii) in the Stokes boundary layer close to the wall. The last contribution has been shown to dominate the other ones 3 . Denoting by µ and ν the water dynamic and kinematic viscosities, and by Ẋ the oscillation speed, we model the viscous dissipation as a single force F w ∼ 2πRλµ Ẋ/δ, where δ ∼ ν/ω 0 is the Stokes layer thickness 3 .

FIG. 4 .

 4 FIG. 4. (A)Numerical nonlinear solution of eq. (3) for h = 4 mm (black line), together with its envelope obtained by asymptotic analysis (orange line). The vertical red line indicates the finitetime of arrest t * . The blue line represents the envelope of the linear model solution starting from the same initial condition. (B) Time of arrest as a function of the foam height h for a single initial condition (from eq.(6) with α(h) reported in Fig.2B). Symbols correspond to the values deduced from the experiments.

Financial support by ERC grant SimCoMiCs 280117 is gratefully acknowledged.

The authors warmly thank Dr Gabs brewery for their inspiring foam, and for providing the fuel for our experiments.


Natacha Ruchaud
email: ruchaud@eurecom.fr

Grigory Antipov
email: grigory.antipov@orange.com

Pavel Korshunov
email: pavel.korshunov@idiap.ch

Jean-Luc Dugelay

Touradj Ebrahimi
email: touradj.ebrahimi@epfl.ch

Sid-Ahmed Berrani

The impact of privacy protection filters on gender recognition

Keywords: Gender Recognition, Deep Learning, Crowdsourcing

Deep learning-based algorithms have become increasingly efficient in recognition and detection tasks, especially when they are trained on large-scale datasets. Such recent success has led to a speculation that deep learning methods are comparable to or even outperform human visual system in its ability to detect and recognize objects and their features. In this paper, we focus on the specific task of gender recognition in images when they have been processed by privacy protection filters (e.g., blurring, masking, and pixelization) applied at different strengths. Assuming a privacy protection scenario, we compare the performance of state of the art deep learning algorithms with a subjective evaluation obtained via crowdsourcing to understand how privacy protection filters affect both machine and human vision.

INTRODUCTION

Recent years show a steady growth in the adoption of digital video surveillance systems for monitoring buildings and public spaces. These systems raise many concerns related to the privacy rights of the subjects being monitored. [START_REF] Dufaux | A framework for the validation of privacy protection solutions in video surveillance[END_REF][START_REF] Fradi | Contextualized privacy filters in video surveillance using crowd density maps[END_REF][START_REF] Moon | Implementation of the privacy protection in video surveillance system[END_REF] At the same time, video analytic tools have created additional challenges, as they can extract privacy sensitive information such as identity, ethnicity, and gender. For instance, face recognition or person reidentification can potentially expose the identity of any individual that appears in the field of view of a wearable camera.

In order to protect privacy of people, several filters have been designed and used, e.g., blurring, masking, or pixelization. The goal of privacy protection filters is to reduce detection and recognition rate of privacy sensitive information. Therefore, privacy filters should prevent the recognition of personal information by humans as well as by computer vision algorithms.

In this paper, we focus on the task of gender recognition in images and evaluate via subjective and objective experiments how the gender recognition abilities of both human vision and computer vision are affected by different privacy protection filters. Previously, humans and computers were compared in their abilities to recognize objects. [START_REF] Borji | Human vs. computer in scene and object recognition[END_REF] It was shown that humans can easily recognize the gender (above 95% of recognition accuracy from faces [START_REF] Bruce | Sex discrimination: how do we tell the difference between male and female faces?[END_REF] ) even in the presence of noise, but such task remains challenging for computer vision. Based on the common knowledge of the world, humans can guess the gender even if an image has been degraded. However, a computer vision model is usually trained only on high-quality images, so the answer remains unclear for the algorithms.

In this paper, we perform subjective experiments using crowdsourcing, which is a viable alternative to conventional laboratory-based subjective assessments. [START_REF] Korshunov | Crowdsourcing-based evaluation of privacy in hdr images[END_REF] Crowdsourcing also provides access to large groups of subjects (or workers) with varying social and geographical backgrounds for a relatively low cost. To display images of body silhouettes to different workers and to collect evaluation results, QualityCrowd2 * framework was used 7 together with Microworkers † crowdsourcing platform, which provided online workers from around the world.

Objective evaluation is done using Convolutional Neural Networks (CNN), which has recently been employed with success on large-scale datasets. We employ Histogram of Oriented Gradients (HOG) features, as they show to be suitable for pedestrian gender recognition. [START_REF] Cao | Gender recognition from body[END_REF][START_REF] Bourdev | Describing people: A poselet-based approach to attribute classification[END_REF] For instance, 76% recognition rate is claimed on MIT dataset, 8 while in Ref. 10, authors announce 80% of the recognition accuracy on the VIPeR dataset.

The images for both crowdsourcing and objective experiments were selected from PETA collection of datasets, 11 which contains 19'000 annotated pedestrian images from various surveillance cameras. Since, in video surveillance, it is difficult to get a clear face of individuals, often due to low quality capture or occlusions, we estimate the gender from a global view of the body silhouette. Methods such as Ellipse fitting [START_REF] Martin-Felez | Towards a more realistic appearance-based gait representation for gender recognition[END_REF] and Radon Transform Mean Gait Energy Image plus Zernike moments [START_REF] Oskuie | Gender classification using a novel gait template: radon transform of mean gait energy image[END_REF] have been implemented for gender recognition from gait.

The rest of the paper is organized as follows. The previous studies related to our work are introduced in Section 2. The datasets and their modifications are presented in Section 3. Then the used gender recognition approaches are described in details and compared in Section 4. Finally, we conclude and give an outlook on possible future works.

BACKGROUND AND RELATED WORK

In this section, we describe privacy protection filters used to generate data for objective and subjective evaluations, as well as, overview the state of the art on gender recognition and crowdsourcing approaches commonly used for subjective assessment.

Privacy filters

To protect privacy, some surveillance systems apply traditional techniques such as masking, Gaussian blurring, and pixelization (e.g., Google Street View, FacePixelizer ‡ , and ObscuraCam § on Android).

More advanced privacy filters have been proposed to hide identity, such as morphing [START_REF] Korshunov | Towards optimal distortion-based visual privacy filters[END_REF] which merges two images and K-means [START_REF] Fradi | Privacy protection filter using shape and color cues[END_REF] which groups colors inside a picture.

Every privacy filter contains a varying parameter which controls the level of privacy protection. The experiments in Ref. 16 demonstrate that, in general, an increase in strength of the filters leads to an increase in privacy (i.e., reduction in recognition accuracy) and at the same time results in a reduced intelligibility (i.e., reduction in detection accuracy).

Masking filter is obtained by applying the following formula:

ImgM asking = originalImg * (1 -α), (1)
with α representing the opacity between 0.5 and 0.9 (the bigger α, the stronger is the filter).

The idea behind the morphing filter 17 is to find an average face image between the source and the target faces according to a given interpolation level. The source face corresponds to the face of the individual whose identity must be preserved. The target face is any generic human face. The method divides both images into Delaunay triangles [START_REF] Benson | Morph transformation of the facial image[END_REF] and transforms the vertices of the source image to the vertices of the target image. The pixel intensities are also interpolated with respect to a second parameter. Morphing is compression and format independent. It is also reversible. Its security can be ensured by encrypting the key points and randomizing the interpolation level and the pixel interpolation values for each triangle. However, as the algorithm begins with triangulating the face images, it may fail to work in cases where the faces are not captured from 'ideal' angles.

Pixelization can be perceived as a downsampling of the images, reducing the overall image quality to protect privacy. The strength of pixelization is controlled by the size of the squares, i.e., ensembles of adjacent pixels (see Figure 1).

Gaussian blurring uses Gaussian kernel to filter an original image:

G(x, y) = 1 2πσ 2 e -x 2 +y 2 2σ 2 , (2)
where x and y are the pixel coordinates and σ is the standard deviation of the Gaussian distribution. The goal of the paper is to estimate the robustness of a gender recognition algorithm against privacy filters and to provide a subjective evaluation based on crowdsourcing. This will also allow us to compare the impact of privacy filters on machine vision and on human vision respectively. The following five filters are considered: black masking with opacities equal to 0.5, 0.7, and 0.9 using Eq. 1, morphing with opacities equal to 0.4, 0.7, and 0.9, pixelization with size of the squares equal to 3, 5, and 7, Gaussian blurring with standard deviations equal to 2, 4, and 6, and K-means with number of clusters equal to 2, 4, and 6. An illustration of these filters is shown in Figure 1 and the used strengths are summed up in Table 5.

Gender recognition in computer vision

Automatic gender recognition might be useful in many applications such as human-computer interaction systems, surveillance systems, content-based indexing and search, biometrics, demographic collection and targeted advertising . [START_REF] Ng | Vision-based human gender recognition: A survey[END_REF] The characteristics of a person (age, ethnicity, facial expressions or accessories) and the image capture process (pose, lighting, illumination, image quality) affect the ability of a computer vision system to recognize the gender.

Guo et al. [START_REF] Guo | Is gender recognition affected by age?[END_REF] report that gender classification accuracy is significantly affected by age (adult faces are recognized more easily than young or senior faces). In our paper, we demonstrate that image quality controlled by a filter (blurring, pixelization, masking, K-means and morphing) can affect the recognition by humans as well as by automatic recognition systems.

Even if our work focuses only on gender recognition based on body silhouettes, in this section, we provide a short overview of gender recognition from other modalities as well. [START_REF] Xia | Multi-view gender classification based on local gabor binary mapping pattern and support vector machines[END_REF] LGBMP hist SVM-RBF 94.96 Wang, 2010 [START_REF] Wang | Dense sift and gabor descriptors-based face representation with applications to gender recognition[END_REF] SIFT, Gabor Adaboost 97 Wu, 2011 [START_REF] Wu | Multi-view gender classification using symmetry of facial images[END_REF] LGBP SVM-RBF 97 Zheng, 2011 [START_REF] Zheng | A support vector machine classifier with automatic confidence and its application to gender classification[END_REF] LGBP-LDA SVMACe 99.8

Table 1. Gender recognition from face.

Gender recognition from face

Table 1 summarizes a list of representative works on face gender recognition.

In many approaches, an image preprocessing (e.g. contrast and brightness normalization or geometric alignment) is applied after a face is detected from an image. Feature extraction methods for face gender classification can be categorized into geometric-based and appearance-based methods. [START_REF] Benabdelkader | A local region-based approach to gender classi. cation from face images[END_REF][START_REF] Li | Gender classification by combining clothing, hair and facial component classifiers[END_REF] Zheng et al. [START_REF] Zheng | A support vector machine classifier with automatic confidence and its application to gender classification[END_REF] obtain the best result for face images from FERET dataset. However, only frontal faces from the dataset were used.

Gender recognition from gait

Works, that are listed in Table 2, are focused on gender recognition from gait (i.e walking, running, jogging and climbing stairs).

First Author, Year Feature extraction

Classifier Average accuracy % Chen, 2009 [START_REF] Chen | Gender recognition from gait using radon transform and relevant component analysis[END_REF] Radon transform + Relevant component analysis Mahalanobis dist 95.7 Chang, 2010 [START_REF] Chang | Using gait information for gender recognition[END_REF] DCT EHMM 94 Felez, 2010 [START_REF] Martin-Felez | Towards a more realistic appearance-based gait representation for gender recognition[END_REF] Ellipse fittings SVM-linear 94.7 Hu, 2010 [START_REF] Hu | Combining spatial and temporal information for gait based gender classification[END_REF] Gabor + MMI GMM-HMM 96.77 Oskuie, 2011 [START_REF] Oskuie | Gender classification using a novel gait template: radon transform of mean gait energy image[END_REF] RTMGEI + Zernike moments SVM 98.9 Oskuie et al. [START_REF] Oskuie | Gender classification using a novel gait template: radon transform of mean gait energy image[END_REF] report the best results for the CASIA gait dataset.

Gender recognition from body silhouettes

One of the first attempts to recognize gender from images of body silhouettes was performed by Cao et al. [START_REF] Cao | Gender recognition from body[END_REF] who used HOG features with Adaboost classifier. Then Collins et al. [START_REF] Collins | Full body image feature representations for gender profiling[END_REF] proposed their descriptor called PixelHOG using dense HOG features computed from an edge map. Biologically-inspired features (BIF) were also used for human body gender recognition by Guo et al. [START_REF] Guo | Gender from body: A biologically-inspired approach with manifold learning[END_REF] Existing methods on gender recognition from body are summarized in Table 3. In Ref. 36 authors train a CNN for gender recognition from body silhouettes on the MIT pedestrian dataset having only 900 images and reach 80 % recognition rate. In Ref. 37 authors obtain the state-of-the-art performance using a CNN trained on PETA collection of datasets.

First

Convolutional neural networks (CNN)

Basic notions in visual neuroscience [START_REF] Hubel | Receptive fields, binocular interaction and functional architecture in the cat's visual cortex[END_REF] led to the creation of the convolutional and pooling layers in CNN. From 1990s, CNN are used for object detection in natural images, including faces and hands [START_REF] Vaillant | Original approach for the localisation of objects in images[END_REF] and for face recognition. [START_REF] Lawrence | Face recognition: A convolutional neural-network approach[END_REF] In images, deep neural networks exploit the property of local combinations of edge motifs which can be regrouped into parts. Firstly, local groups of pixels are often highly correlated in images, forming distinctive local motifs that are easily detected. Secondly, the local statistics of images are invariant to location.

The architecture of a typical CNN is structured as follows: the first few stages are composed of two types of layers: convolutional layers and pooling layers. Units in a convolutional layer are organized in feature maps and each unit is connected to local patches in the feature maps of the previous layer through a set of weights called a filter bank. The result of this filter bank is then passed through a non-linearity such as a ReLU. The same filter bank is connected with all units in a feature map. Mathematically, a discrete convolution is the filtering operation performed by a feature map.

The role of the convolutional layer is to detect local correlations of features from the previous layer and the role of the pooling layer is to merge semantically similar features into one. The maximum of a local patch of units in one feature map is computed in a typical pooling unit. Usually two or three stages of convolution, non-linearity and pooling are cumulated. Finally, back-propagation gradients allow all the weights in all the filter banks to be trained.

Privacy assessment based on crowdsourcing

Until recently, little was done to better understand privacy issues in practical multimedia applications. But lately the impact of privacy protection tools has been analyzed in video surveillance and effective evaluation methodologies have been developed to take into account both the context and the content. The objective evaluation of several primitive privacy filters has been reported by Newton et al., [START_REF] Newton | Preserving privacy by de-identifying face images[END_REF] where the authors showed that such filters cannot adequately protect from successful face recognition, because recognition algorithms are robust. The robustness of face recognition and detection algorithms to primitive distortions is also reported in Ref. 42. Further, in a work by Dufaux et al., [START_REF] Dufaux | A framework for the validation of privacy protection solutions in video surveillance[END_REF] a framework is defined to evaluate the performance of face recognition algorithms applied to images altered by various obfuscation methods.

Crowdsourcing has shown to be a viable alternative to conventional laboratory-based subjective assessments, especially for cognitive tasks. [START_REF] Korshunov | Crowdsourcing-based evaluation of privacy in hdr images[END_REF] Crowdsourcing-based evaluation of privacy tradeoff in video surveillance has shown good consistency with laboratory-based studies. [START_REF] Korshunov | Crowdsourcing approach for evaluation of privacy filters in video surveillance[END_REF] The crowdsourcing methodology benefits from a large number of participants and can be performed efficiently and at a relatively low cost without requiring a significant commitment from subjects, which are called workers in the crowdsourcing terminology. Workers accept to undertake a task (usually a short 5-20 minutes task) and are grouped in larger units, called batches. When the evaluation experiment is over, workers submit their answers. Unlike laboratory-based experiments, crowdsourcing cannot impose specific displays or controlled illumination of surroundings in which assessments take place. However, since standard environment and equipment conditions for surveillance operators have not been established, typical monitors even with different resolutions and color settings are considered as appropriate in this study.

To display video sequences to different workers and to collect evaluation results, we selected QualityCrowd2 framework [START_REF] Keimel | Qualitycrowd a framework for crowd-based quality evaluation[END_REF] and the Microworkers crowdsourcing platform in order to access online workers from around the world. QualityCrowd2 is an open-source framework designed for QoE evaluation with crowdsourcing. This framework was selected because it is easy to modify for our gender recognition task using the provided simple scripting language for batch creation and training sessions.

DATASET DESCRIPTION

Contrary to previous studies, [START_REF] Cao | Gender recognition from body[END_REF][START_REF] Collins | Full body image feature representations for gender profiling[END_REF] we do not focus on a single dataset, but on a collection of pedestrians datasets PETA. [START_REF] Deng | Pedestrian attribute recognition at far distance[END_REF] An extensive evaluation of a recognition algorithm requires a dataset containing various people in diverse environments and situations. Therefore, a large and diverse database of pedestrians has been selected. PETA is a large open-access collection of pedestrian images with several annotations including age, backpack, hat, jacket, jeans, logo, long hair, gender, muffler, no accessory, plastic bag, sandals, shorts, skirt, sunglasses, trousers, t-shirt, etc. Examples of PETA images are presented in Figure 2. Originally, the PETA collection consists of 10 datasets of different sizes with a total of 19'000 images. Appearances of images significantly vary between different datasets of PETA in terms of image resolutions (from 17x39 to 169x365 pixels), camera angles (pictures are taken either by ground-based cameras or by cameras which are set at a certain height), and environments (indoors or outdoors), as shown in Figure 3. Thanks to the proposed annotations, we can evaluate gender recognition accuracy on the PETA collection of datasets.

Unfortunately, PETA contains many images of the same persons which can considerably bias the resulting prediction rate. This is why these repetitions are removed from PETA. Moreover, images with very low resolutions (when height is less than 120 pixels or width is less than 40 pixels), images with more than one person, and images of babies in strollers are also removed. Finally, we are left with 8'365 images which is less than half of the initial size of PETA.

EVALUATION EXPERIMENTS 4.1 CNN-based gender recognition

Following the work in Ref. 37, we employ a convolutional neural network for gender recognition. In particular, we adopt an architecture proposed by Krizhevsky et al. [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF] This architecture is presented in Figure 4. It consists of five convolutional layers and three fully-connected layers. We train a variant of this architecture on PETA dataset. The variation is in the last fully-connected layer, where we use two neurones instead of 1'000, since we only have two target classes (male and female). The network was trained using the Tesla K20c graphical processor. In order to generalize our CNN on heterogeneous and even completely unseen datasets, we firstly train CNN on the dataset which is composed of training parts of CUHK, PRID, GRID, MIT, and VIPeR taken together (4488 Male and 2744 Female pedestrian images). Then, our model is tested on a single big dataset composed of images from the 10 datasets which have not been used in the training set (160 Male and 161 Female pedestrian images). Table 4 summarizes the exact number of training and testing images for each dataset.

Privacy protection filters at different strengths, as described in Table 5, have been applied only on the testing part.

Dataset

Train size (male + female) Test size (male + female) CUHK 4. Split between training and testing parts per dataset.

In order to evaluate performance of the model, the mean average precision (MAP) is used. Figure 5 shows that the trained CNN is robust to K-means and pixelization (more than 70% of recognition accuracy even if the filter impact is strong), while it is sensible to morphing (less than 60% of recognition accuracy).

Crowdsourcing Evaluation

The crowdsourcing assessment aims at checking if a person in a given image can be correctly identified as female or male by an individual, even after the privacy protection For this purpose, each crowdsourcing worker was asked to look at the image of a person and answer the question "What is the gender of the person?" with the following options: "male", "female", and "I don't know".

In total, 300 random images from the PETA dataset were used in the crowdsourcing experiment. They were protected by five different privacy filters at three different strength levels, resulting in 300 × 5 × 3 = 4500 images evaluated in this crowdsourcing study.

To ensure a statistically significant number of evaluations for each image, also taking into account the presence of unreliable subjects (about 50% in a typical crowdsourcing evaluation), 40 subjects were assigned to each image, with a total of 2652 subjects participating in the evaluations.

All versions of the images were randomly distributed among the batches; special care was devoted to guarantee that a particular content was used only once in each batch, i.e., that each subject assessed only one version of a given content. Each batch starts with a message about what is required from the subject, followed by a training session describing the evaluation procedure. A display brightness test is performed using a method similar to that described in Ref. [START_REF] Hossfeld | Best practices for QoE crowdtesting: QoE assessment with crowdsourcing[END_REF] and permits to estimate the subjects' display settings. Subjects are not allowed to skip any content or to avoid answering any question.

Unlike laboratory-based subjective experiments where all subjects can be observed by operators and its test environment also can be controlled, the major shortcoming of the crowdsourcing-based subjective experiments is the inability to supervise participants behavior and to restrict their test conditions. When using crowdsourcing for evaluation, there is a risk of including untrusted data into analysis due to the wrong test conditions or unreliable behavior of some workers who try to submit low quality work in order to reduce their effort while maximizing their received payment. [START_REF] Hossfeld | Best practices for QoE crowdtesting: QoE assessment with crowdsourcing[END_REF] For this reason, unreliable workers detection is an inevitable process in crowdsourcing-based subjective experiments. To identify a worker as 'trustworthy', the following four factors were taken into account in our experiments:

• Task completion time;

• Mean observation time per question;

• Observation duration deviation;

• Number of minority answers.

The objective of the first three factors is to filter out the workers who have strange behaviors in the middle of their tasks, because they are either not serious or have poor concentration. The observation time per question is measured as the time from when the question is displayed until the time the answer is given by the worker. The task completion time, mean observation time and observation duration deviation can be calculated using this data. If the task completion time or mean observation time per question is too long compared to the average of all workers, it can be deduced that they did not take the test seriously or were distracted during their tasks. Unreliable workers were also identified using an approach similar to a typical outlier detection method, commonly used in most subjective quality evaluations. However, typical subjective tests use scoring methods like five-grade evaluation, and outlier detection is performed on mean opinion score. [START_REF]Methodology for the subjective assessment of the quality of television pictures[END_REF] Our experiments do not have opinion scores because of the specific privacy-oriented questions we used. Therefore the number of minority answers has been used for the outlier detection instead. The assumption is that a participant who has a lot of different answers when compared to the majority of workers is unreliable. The threshold was determined in the same way as for the other three factors.

The above unreliable worker detection methods filtered out 198 workers out of total 2652, resulting in 2454 scores used in the analysis.

Results

An output of the CNN is binary ("male" or "female") while the output of the crowdsourcing is ternary ("male", "female" "I don't know"). Therefore, in order fairly compare the results of the gender recognition based on CNN and crowdsourcing, we assume that 50% of answers "I don't know" registered during the crowdsourcing had been correct if the response would have been selected randomly. Figure 5 illustrates the results of the CNNbased gender recognition as well as crowdsourcing results. Results are given for filtered images (with different parameter values) and for original images (i.e. no filter). The CNN shows results close to those by crowdsourcing for Gaussian blurring, K-means and morphing. CNN is more robust to pixelization (∼ 10% better than human vision). In the case of masking, the human vision is ∼ 15% better than the CNN. Protection of the masking filter really depends on the brightness of the display and its environment. Indeed, some crowdsourcing participants might have increased the luminosity of their computer screens making the test images more visible, whereas the CNN works with values of pixels regardless of the display method.

CONCLUSION AND FUTURE WORK

This paper demonstrates the impact of privacy filters on gender recognition by machine vision algorithms and by human vision. We assessed gender objectively (using a CNN) and subjectively (using a crowdsourcing approach). One might expect human vision to be more robust to privacy filters than computer vision. Nevertheless, our results show that humans and automatic gender recognition systems perform almost equally.

In our future work, we would like to explore and compare the nature and causes of errors committed by humans and automatic gender recognition systems.

  The K-means filter works with pixel values. K-means minimizes the Euclidian distance between a pixel and a centroid of the closest color. First, a number of clusters (denoted here by N) is chosen. Then K-means clustering is applied on the color space (r, g, b) and N centroids are obtained. Finally, original pixel values are replaced by the corresponding centroid values. Strength of K-means is controlled by the number of clusters.

Figure 1 .

 1 Figure 1. Privacy filters. From left to right: original body, masking, morphing, pixelization, Gaussian blur, K-means.

Figure 2 .

 2 Figure 2. PETA dataset.

Figure 3 .

 3 Figure 3. PETA dataset.

Figure 4 .

 4 Figure 4. CNN model.

Figure 5 .

 5 Figure 5. Accuracy results of human vision and CNN.

Table 2 .

 2 Gender recognition from gait.

Table 3 .

 3 Gender recognition from body silhouettes.

Table 5 .

 5 Privacy filters and the corresponding varying parameters.

	Filter	Parameter	Parameter values
	Black Masking opacity	0.5, 0.7, 0.9
	Morphing	opacity	0.4, 0.7, 0.9
	Pixelization	size of squares	3, 5, 7
	Gaussian Blur	standard deviation 2, 4, 6
	Kmeans	number of cluster	6, 4, 2

ACKNOWLEDGEMENTS

This work has been partially performed under the European Network of Excellence VideoSense under FP7 ¶ . ¶ http://videosense.eu/