David A Bonhommeau
email: david.bonhommeau@univ-reims.fr

MCMC 2 (version 1.1.1): A Monte Carlo code for multiply charged clusters

Keywords: .gz. Programming language: fortran 90 with MPI extensions for parallelization. Computers: x86 and IBM platforms Monte Carlo simulations, Coarse-grained models, Charged clusters, Charged droplets, Electrospray ionisation, Parallel Tempering, Parallel Charging

This new version of the MCMC 2 program for modeling the thermodynamic and structural properties of multiply-charged clusters fixes some minor bugs present in earlier versions. A figure representing the required RAM per replica as a function of the cluster size (N ≤ 20000) is also provided as benchmark.

Has the code been vectorised or parallelized? : Yes, parallelized using MPI extensions. Number of CPUs used: up to 999. RAM (per CPU core): 10-20 MB. The physical memory needed for the simulation depends on the cluster size, the values indicated are typical for small or medium-sized clusters (N ≤ 300 -400). The size of A n+ N clusters (N = number of particles, n = number of charged particles with n ≤ N ) should not particular, in test run 2, the particle and probe diameters used to evaluate the number of surface particles were fixed to 0.8 and 1.2, respectively (see keyword "SURFACE"). Actually, the probe diameter should be smaller than the particle diameter [START_REF] Miller | Structure and stability of charged clusters[END_REF] and the two values were therefore swapped. (c) The subroutines dLJ nopol hom (in lib4-nopol.f90), dLJ pol hom (in lib4-pol.f90), and dLJ dampol hom (in lib4-dampol.f90) are renamed dLJ nopol, dLJ pol, and dLJ dampol, respectively, to avoid any ambiguity. The suffix "Hom", that stood for "homogeneity" in order to indicate that Lennard-Jones interactions between particles were the same, was improper since homogeneity is commonly related to invariance by translation and all the properties of multiply charged clusters cannot be considered invariant by translation in the most general case. The renaming of the three subroutines has obviously no influence on the results and some related comments have been modified accordingly.

Acknowledgements

Dr Mark A. Miller and Prof. Marius Lewerenz are gratefully acknowledged for providing some routines and valuable advices during the development of the program. This work was supported by the "Centre National de la Recherche Scientifique" (CNRS) through an excellence chair (DAB).

exceed 1.6 × 10 4 (respectively 2.0 × 10 4 ) particles on servers with 2 GB (respectively 3 GB) of RAM per CPU core if n = 0 (neutral clusters) or n = N ("fully-charged" clusters). For charged clusters composed of neutral and charged particles (eg, n = N/2), the maximum cluster size can drop to 1.4 × 10 4 and 1.8 × 10 4 on servers with 2 GB and 3 GB of RAM, respectively (see the figure given in Supplementary Material).

Nature of problem:

We provide a general parallel code to investigate structural and thermodynamic properties of multiply charged clusters.

Solution method : Same as in the previous version.

Restrictions:

The current version of the code uses Lennard-Jones interactions, as the main cohesive interaction between spherical particles, and electrostatic interactions (charge-charge, charge-induced dipole, induced dipole-induced dipole, polarisation). Furthermore, the Monte Carlo simulations can only be performed in the N V T ensemble and the size of charged clusters should not exceed 2.0×10 4 particles on CPU cores with less than 3GB of RAM each. It is worth noting that the latter restriction is not significantly crippling since MCMC 2 should be mainly devoted to the investigation of medium-sized cluster properties due to the difficulty to converge Monte Carlo simulations on large systems (N ≥ 10 3 ) [START_REF] Pahl | Accurate Melting Temperatures for Neon and Argon from Ab Initio Monte Carlo Simulations[END_REF].

Unusual features: Same as in the previous version.

Running time: Same as in the previous version.

Summary of revisions

1. Additional features of MCMC 2 version 1.1.1: Same as in the previous version.

2. Modifications or corrections to MCMC 2 version 1.1 [START_REF] Bonhommeau | MCMC 2 : A Monte Carlo code for multiply-charged clusters[END_REF][START_REF] Bonhommeau | MCMC 2 (version 1.1): A Monte Carlo code for multiply-charged clusters[END_REF] (a) Several minor bugs were fixed in this version i. A default value for the integer "irand", used to select the type of random number generator (keyword SEED, subkeyword METHOD), was missing. It is set to 0. ii. The subkeyword "EVERY" used to define the frequency of statistics printing (keyword "STATISTICS") was missing and it has been implemented in the program. Before version 1.1.1, the choice entered into the setup file was simply ignored and the frequency was always set to its default value, namely a printing every 100 Monte Carlo sweeps. (b) Some useless integers are removed from subroutines in lib4-pol.f90 and lib4-dampol.f90 and some test runs are slightly modified. In


David A Bonhommeau
email: david.bonhommeau@univ-reims.fr

MCMC 2 (version 1.1.2): A Monte Carlo code for multiply charged clusters

Keywords: .gz. Programming language: fortran 90 with MPI extensions for parallelization. Computers: x86 and IBM platforms Monte Carlo simulations, Coarse-grained models, Charged clusters and droplets, Electrospray ionisation

This version of the MCMC 2 program corrects several glitches of earlier versions and adapts the keyword environment of radial, angular, and evaporation histograms to remove some useless columns in the data files or simply prevent their creation, when desired. Although most of the fixed bugs only resulted in warnings at compiling, one of them, related to an incorrect array reinitialization, was identified to be the main source of instabilities of earlier versions of the program when compiled with the ifort compiler. The setup keywords altered to take into account the aforementioned corrections and improvements are detailed at the end of the manuscript.

Has the code been vectorised or parallelized? : Yes, parallelized using MPI extensions. Maximum number of CPUs allowed: 999.

RAM (per CPU core): 10-20 MB. The physical memory needed for the simulation depends on the cluster size, the values indicated are typical for small or medium-sized clusters (N ≤ 300 -400). The limitations on the neutral and charged cluster sizes are the same as in MCMC 2 version 1.1.1 [START_REF] Bonhommeau | MCMC 2 (version 1.1.1): A Monte Carlo code for multiply charged clusters[END_REF].

Nature of problem:

We provide a general parallel code to investigate structural and thermodynamic properties of multiply charged clusters.

Solution method : Same as in the previous version.

Restrictions: Same as in the previous version.

Unusual features: Same as in the previous version.

Running time: Same as in the previous version.

Reasons for the new version:

Several glitches that may result in warnings at compiling or instabilities at the execution stage have been corrected. In particular, one of these glitches related to an incorrect array reinitialization might lead to crashes of the program. This bug is deemed to be the source of the instabilities denoted in earlier versions of the code when they were compiled with the ifort compiler. We also take advantage of this new version to improve the definition of the keywords devoted to the generation of histograms and we correct some comments throughout the program files. i. The arrays racc and rrej used to evaluate the acceptance and rejection rates on each replica after PT or PC configuration swap attempts are converted into double precision reals since only one component of these arrays needs to be known per replica. This should save a small amount of memory. Amended subroutines: simulate nopol in lib4-nopol.f90, simulate pol in lib4-pol.f90, simulate dampol in lib4-dampol.f90, ptswap, pcswap nopol, pcswap pol, pcswap dampol, ptswap2, pcswap2 nopol, pcswap2 pol, and pcswap2 dampol in lib5.f90. ii. A comment in the main file (main.f90) was incomplete: The second stage of the equilibration phase can be PT or PC depending on the setup file entries, and not PT only. iii. The dimension of revap, array that stores the number of evaporated particles, was incorrectly fixed to 3 instead of 5 in printinstats nopol (lib3-nopol.f90) which altered the intermediate files created at the end of each stage of a split simulation [START_REF] Bonhommeau | MCMC 2 : A Monte Carlo code for multiply-charged clusters[END_REF] involving charged particles. iv. The dimension of N pol and N nopol arrays is set to 4 instead of 2 in printout evap pol (lib3-pol.f90) to comply with their initial declaration in the module config (lib2.f90) although the two last components of these vectors are not currently used in printout evap pol. (d) Improvement of the evaporation histograms.

Summary of revisions

Additional features of MCMC

i. A logical levap is added to enable the user not to generate the evaporation histograms. ii. A logical levap2 is added to enable the user to reinitialize to zero the evaporation histograms during split simulations. iii. Some comments are now printed in the output file to indicate whether the generation of evaporation histograms is deactivated. Amended subroutines and modules: run params in lib2.f90, read setup in read setup.f90, evap nopol in lib3-nopol.f90, evap pol in lib3-pol.f90, simulate nopol in lib4-nopol.f90, simulate pol in lib4-pol.f90, and simulate dampol in lib4-dampol.f90. (e) Improvement of the angular distributions. The angular distributions represent the angles between the radius vectors of charge carriers. In earlier versions of the code the angular distributions were systematically evaluated for all the charge carriers within the cluster and for the charge carriers located at the cluster surface, whatever the cluster composition. In particular, the angular distribution files were filled with zeroes when the cluster was neutral. In this case, the calculation of the angular distributions is now automatically deactivated by setting the input parameter cos every to zero. A warning is added in the output file to inform the user of this reinitialization. Amended routines and subroutines: main program MCMC2 in main.f90, print surf, printin surf and readin surf in lib3.f90, histo surf nopol in lib3-nopol.f90, histo surf pol in lib3-pol.f90, simulate nopol in lib4nopol.f90, simulate pol in lib4-pol.f90, and simulate dampol in lib4dampol.f90. (f) Improvement of the radial histograms.

i. We remove the possibility to set radtyp = 0, where radtyp is the parameter to select the type of radial histograms to be computed. This parameterization was redundant with the subkeyword USE of keyword RADIAL that was introduced in version 1.1 [START_REF] Bonhommeau | MCMC 2 (version 1.1): A Monte Carlo code for multiply-charged clusters[END_REF] to enable the deactivation of the calculation of radial histograms.

The remaining integral values allowed for radtyp, namely 1 -4, are subsequently shifted to 0 -3. ii. radtyp = 1 (ie, generation of radial distributions for charged particles) is forbidden when the cluster is neutral. radtyp is then reinitialized to zero (default value). iii. radtyp = 2 (ie, generation of radial distributions for neutral particles) is forbidden when all the particles are charged. radtyp is then reinitialized to zero. iv. radtyp = 3 (ie, generation of radial distributions for all the particles, for charged particles, and for neutral particles) is forbidden when the cluster is neutral, only composed of positively charged particles, or only composed of negatively charged particles. In these three cases radtyp is reinitialized to zero. These reinitializations of radtyp aimed to prevent the production of useless columns in the data files. Their omission in previous versions of the MCMC 2 program did not affect the accuracy of the MC results. Amended routines and subroutines: main program MCMC2 in main.f90, read setup in read setup.f90, inhisto, inhisto split, normhisto, printin histo, readin histo in lib3.f90, calhisto nopol in lib4-nopol.f90, calhisto pol in lib4-pol.f90 (g) Modified comments. Some comments and column headers of the data files have been corrected.

Modified setup keywords

• RADIAL USE lrad GRIDRCOM deltagrid stepgrid PARTICLE radtyp: Setting parameters for generating one-particle radial histograms. lrad indicates whether the user wants (lrad = .true.) or does not want (lrad = .false.) the calculation of radial histograms. These histograms cannot extend beyond the radius of the Monte Carlo container but, for graphical purposes, the user is allowed to add a small distance deltagrid to the grid size (= radius + deltagrid). stepgrid is the grid step for one-particle radial histograms with respect to the cluster center of mass. The grid origin is hardcoded to zero since these histograms are calculated with respect to the cluster center of mass. The number N grid of grid points is automatically determined in the code from the grid size and the grid step. radtyp enables the user to specify the type of particles to be considered in the radial histograms (0 = all the particles without any distinction, 1 = charged particles only, 2 = neutral particles only, 3 = all the previous histograms (3=0+1+2)). In case of inconsistency between the input configurations and the value of radtyp, this parameter is reinitialized to zero. Default: lrad = .false., deltagrid = 0, stepgrid = 0.1, radtyp = 0.

• EVAPORATION USE levap PARTNB N part evap GRIDNB N samp evap: Setting of parameters for defining the evaporation criterion and related histograms. levap indicates whether the user wants (levap = .true.) or does not want (levap = .false.) the calculation of evaporation histograms. N part evap defines the number of particles above which evaporation is assumed to be prominent during a MC sweep. N samp evap is the bin size (in number of particles) selected for plotting the distribution of evaporated particles. Default: levap = .false., N part evap = 1 and N samp evap = 1.

• SPLITSIM USE lsplit isplit nsplit ANALYSIS lsurf lrad2 lang levap2: Setting of parameters for running split MC simulations. lsplit indicates whether the user wants (lsplit = .true.) or does not want (lsplit = .false.) to run a split simulation. Each split simulation is composed of several stages where isplit is the number of the current stage and nsplit is the total number of stages. lsurf , lrad2, lang, and levap2 enable the user to reset the arrays for surface particle distributions, radial histograms, angular distributions of charged particles, and evaporation histograms (lsurf = .true., lrad2 = .true., lang = .true., and levap2 = .true.) or to use intermediate files (ie, output files from the previous stage) as starting point of the next stage (lsurf = .false., lrad2 = .false., lang = .false., and levap2 = .false.). Default: lsplit = .false., isplit = 1, nsplit = 2, lsurf = .false., lrad2 = .false., lang = .false., levap2 = .false.

  Correction of minor glitches. Recent gfortran and ifort compilers might generate some warnings at compiling MCMC 2 version 1.1.1. Several subroutines have been slightly modified to remove these warnings. i. Alterations in lib1.f90: The integer iout is no longer defined as a parameter in the subroutine vranf (main routine of the lagged Fibonacci random number generator). It is also declared as INTENT(INOUT) instead of INTENT(IN) in the subroutines xuinit and xuwarm. ii. Alterations in lib4-nopol.f90: A dummy value "0" used in the call of the subroutine dLJ nopol is replaced by an array called dummy (of dimension N part and declared as real) in the subroutine simulate nopol. iii. Alterations in lib4-pol.f90 and lib4-dampol.

	instabilities of the MC simulations attributed to possible issues with
	the ifort compilers [2] were probably a manifestation of this glitch.
	Amended subroutines: simulate nopol in lib4-nopol.f90, simulate pol
	in lib4-pol.f90, and simulate dampol in lib4-dampol.f90.
	(b)
	2. Modifications or corrections to MCMC 2 version 1.1.1:
	(a) Major glitch correction. One major glitch related to an improper
	reinitialization of the logical array lpsim swap devoted to the storage
	of rejection and acceptance decisions after Parallel Tempering (PT)
	or Parallel Charging (PC) configuration swap attempts has been cor-
	rected. In some cases, especially for the replica of higher number,
	the program attempted to reinitialize a component of the array lp-
	sim swap that exceeded the array dimension which might cause mem-
	ory errors. We performed several simulations of 10 10 Monte Carlo
	(MC) steps on small clusters (100 particles) and large clusters (1000
	particles) to ensure that the program was stable when compiling the
	code with several versions of the gfortran and ifort compilers. Former

2 version 1.1.2: Same as in the previous version. f90: A dummy value "0" used in the call of the subroutine dLJ pol (respectively dLJdampol) is replaced by an array called dummy (of dimension N part and declared as real) in the subroutine simulate pol (respectively simulate dampol). The arrays dr (of type vector and dimension N part × N part) and dist (of type real and dimension N part × N part) are declared as INTENT(INOUT) instead of INTENT(IN) in the subroutine LJ pol2 (respectively LJ dampol2) . (c) Other minor corrections.

Acknowledgements

Dr Mark A. Miller, Prof. Marius Lewerenz, and Prof. Marie-Pierre Gaigeot are gratefully acknowledged for providing some routines and valuable advices during the long-term development of the program. DAB is indebted to the computer center of Champagne-Ardenne (ROMEO) for computational resources and technical support.


A Chastagner

M Moinet

Grégoire Perez

E Roy

Karen Mccoy

Olivier Plantard

A Agoulon

S Bastian

Alain Butet

Yann Rantier

Prevalence of Anaplasma phagocytophilum in small rodents in France

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


Romain Casati

Gilles Daviet

Florence Bertails-Descoubes

Inverse Elastic Cloth Design with Contact and Friction

Keywords: CR Categories: I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling -Physically based modeling I.3.7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Animation cloth modeling, frictional contact, inverse problem, nonsmooth optimization

Physically based cloth modeling is classically achieved through a trial and error process. The rest (undeformed) configuration of the cloth, often represented as a 2D pattern assembly, is edited geometrically and adjusted iteratively depending on the feedback provided by a static cloth simulator, which predicts the deformed 3D shape under gravity and contacts. Matching a reference 3D shape while keeping the time of the modeling process reasonable is thus difficult, unless the user possesses advanced skills in real cloth tailoring. In contrast, in this paper we investigate a new, inverse strategy for modeling realistic cloth intuitively. Our goal is to take as input a target (deformed) 3D shape, and to interpret this configuration automatically as a stable equilibrium of a cloth simulator, by retrieving the unknown rest shape. In the presence of gravity and frictional contact, such an inverse problem formulates as an illposed nonlinear system subject to nonsmooth constraints. To select and compute a plausible solution, we design an iterative two-step solving process. In a first step, contacts are reduced to frictionless bilateral constraints, and starting from an as-flat-as possible pose, a unique rest pose is retrieved using the adjoint method on a regularized energy. The second step modifies this rest pose so as to project bilateral forces onto the admissible Coulomb friction cone, for each contact. We show that our method converges well in most cases towards a plausible rest configuration, and demonstrate practical inversion results on various cloth geometries modeled by an artist.

Introduction

Cloth simulation has been a very active research area over the last decades, yielding compelling animations of cloth in movies that are now hardly distinguishable from reality. In parallel, the variety of garments simulated, now ranging from simple cloaks to complex multi-layered dresses, has been considerably enriched in recent years. Such a diversity in shape and motion encourages artists to create personalized garments expressing not only aesthetics and social rank of a given character, but also some more personal traits (e.g., severity, shyness, spirit), which actively contributes to the overall story telling.

However, modeling a specific virtual garment still remains a difficult and tedious task for a designer, as current modeling tools offer Figure 1: Breathing physics into a purely geometric garment directly modeled in 3D under gravity. From the 3D mesh of the cloth and the body mesh, our method automatically infers the rest pose of a cloth simulator such that a stable equilibrium of the cloth simulator under gravity and frictional contact closely matches the input geometry. Starting from this input pose, the cloth can then be animated in a plausible manner, for instance under wind.

very little user control over the final (simulated) shape of the cloth. On the one hand, a number of sophisticated tools allow to tailor garments as one would do in real: The user designs 2D cloth patterns, then sews them together in 3D, and finally adjusts the result iteratively based on the simulation of the "drape" -that is, the stable equilibrium pose of cloth under external forces [START_REF] Volino | From early virtual garment simulation to interactive fashion design[END_REF]; Marvelous Designer 2010; [START_REF] Umetani | Sensitive couture for interactive garment editing and modeling[END_REF]. Though useful to design realistic patterns which may then be manufactured in real, such a trial and error process remains complex and unnatural when dressing a virtual character with a precise final shape in mind (e.g., the designer should anticipate the length of the seams so as to obtain folds of a given size located in a specific area). On the other hand, many different WYSIWYG approaches, ranging from sketch-based modeling [START_REF] Turquin | A sketch-based interface for clothing virtual characters[END_REF] to image-based capture [START_REF] Bradley | Accurate multiview reconstruction using robust binocular stereo and surface meshing[END_REF], have provided the user with intuitive and semi-automatic tools for digitizing the 3D shape of a specific garment, from an imaginary or real reference. Such geometric tools allow the user to focus only on the final shape of the cloth, thus providing a maximal level of control over geometric details. However, although the resulting geometry generally corresponds to an equilibrium configuration of the cloth under external forces, it cannot be interpreted physically in a simple way. If one naïvely plugs such a final pose as the rest (undeformed) state parameter of a cloth simulator, sagging will occur when simulation is activated, and all the modeling efforts and time will be lost.

In this paper we explore a new strategy for modeling virtual cloth both realistically and intuitively, by combining the fine user control offered by geometric tools together with the predictive capabilities of simulation. Our idea is to devise an automatic inverse modeling process able to interpret a merely geometric drape pose as a stable equilibrium of a cloth mechanical simulator. Assuming material properties (mass, stiffness, friction coefficients) of the cloth are known, our algorithm takes as input a target (deformed) 3D surface, resulting from either geometric design or capture, and automatically retrieves a plausible rest shape for the cloth (non-necessarily flat due to seams) as well as frictional contact forces at play. When simulated with such parameters under external forces, the cloth closely matches the input pose, and may be subsequently animated.

Related Work

Cloth modeling and simulation has a rich history in Computer Graphics and Computer-Aided Design. Within less than three decades, practicable solutions have been proposed to design extremely complex garments and to animate them with an impressive level of details. Such remarkable advances were seamlessly integrated into the work flow of the entertainment industry, yielding compelling visual effects in movies and arousing a growing interest in the fashion industry. Here we briefly discuss major approaches for cloth simulation before focusing on current techniques for modeling realistic garments. A recent survey on the topic, including advances in both Computer Graphics and Computer-Aided Design, can be found in [START_REF] Liu | A survey on {CAD} methods in 3d garment design[END_REF].

Cloth Simulation Since the 90's, a large body of work in Computer Graphics contributed to the simulation of cloth dynamics. The work by [START_REF] Baraff | Large steps in cloth simulation[END_REF] was the first to allow for the stable simulation of full-size garments in reasonable timings, thanks to a first-order implicit discretization of the dynamics. Subsequent works contributed on the robust treatment of contact and friction [START_REF] Bridson | Robust treatment of collisions, contact and friction for cloth animation[END_REF], on better formulations for bending [START_REF] Bridson | Simulation of clothing with folds and wrinkles[END_REF][START_REF] Grinspun | Discrete Shells[END_REF]], and on the limitation of stretching [START_REF] Goldenthal | Efficient simulation of inextensible cloth[END_REF][START_REF] English | Animating developable surfaces using nonconforming elements[END_REF][START_REF] Thomaszewski | Continuum-based strain limiting[END_REF]]. Current challenges include the incorporation of real material properties (see, e.g., [START_REF] Wang | Datadriven elastic models for cloth: modeling and measurement[END_REF]) and the question of interactivity at high resolutions (see, e.g., [START_REF] Kim | Near-exhaustive precomputation of secondary cloth effects[END_REF][START_REF] Hahn | Subspace clothing simulation using adaptive bases[END_REF]).

Despite the profusion of simulation models, it is noteworthy that all of them share a similar nodal representation of cloth. That is, cloth is always modeled as a mesh whose vertices (3D positions) serve as degrees of freedom for the corresponding dynamical system.

Physics-based cloth tailoring

The real process for making wearable garments involves the design of 2D patterns, i.e., flat patches which are then manufactured with real fabric and sewn together to create the final garment. Early cloth simulation methods mimicked this process in order to dress virtual characters before animating them [START_REF] Carignan | Dressing animated synthetic actors with complex deformable clothes[END_REF]]. To improve the garment modeling process, Volino and colleagues [2005] have later imagined an intuitive and interactive design environment, where the user can simultaneously edit 2D patterns and visualize the 3D result under gravity and body contacts thanks to a fast draping simulator. Such an interactive physics-based tailoring process is now the standard work flow used in many commercial tools like the popular Marvelous designer software [2010]. It has also inspired further work in academia [START_REF] Umetani | Sensitive couture for interactive garment editing and modeling[END_REF], which has improved the 3D responsiveness to 2D pattern editing and thus greatly smoothed the creative flow.

Virtual tailoring is powerful in the sense that it allows experts in pattern design to create and experiment with complex garments without having to manufacture them. Moreover, final results may not only be used for making real garments, but also for dressing virtual avatars in a realistic way. This has motivated a number of VFX studios to adopt such a modeling strategy [Marvelous Designer 2010;CreativeBloq 2012]. An advanced usage however requires some specific skills in pattern sewing, making the technique pretty complex and unnatural for traditional 3D designers [CreativeBloq 2012]. Moreover, such a trial and error process proves unpractical for reproducing the exact look of a visual reference provided for instance as a painting, a photograph or a 3D capture of a real drape, .

Geometry-based cloth modeling

An interesting alternative in cloth modeling is to leverage 3D geometric tools for directly sculpting the final pose regardless of the underlying physical process at the origin of the cloth deformation. Two families of geometric methods are currently available: manual geometric design, and automatic geometric acquisition. On the one hand, advanced shape editing tools allow 3D designers to carve the drape of a garment directly around a virtual character (see, e.g., [START_REF] Porumbescu | Shell maps[END_REF]). For fairly simple garments, sketch-based interfaces may greatly improve the intuitiveness and speed of the process [START_REF] Igarashi | Clothing manipulation[END_REF][START_REF] Turquin | A sketch-based interface for clothing virtual characters[END_REF]]. Once a final garment has been created, some geometric transfer methods can be used for automatically adapting it to various character morphologies, while preserving the style [START_REF] Brouet | Design preserving garment transfer[END_REF]]. On the other hand, thanks to the considerable advance of image-based capture these latest years, it becomes now affordable to acquire precisely the full 3d geometry of static cloth with folds and wrinkles [START_REF] White | Capturing and animating occluded cloth[END_REF][START_REF] Bradley | Accurate multiview reconstruction using robust binocular stereo and surface meshing[END_REF], thus allowing for an automatic garment creation from a real source.

Geometry-based techniques are appealing because they provide the user with full control over the final 3D shape of the garment. Realistic but also imaginary garments can be created this way, leaving aside the flat patch-based structure of garments and focusing instead on the free-form 3D appearance. Yet, since the modeling process is done independently from physical simulation, the resulting geometry cannot be interpreted mechanically in a straightforward way. This implies that the modeled garment cannot be animated nor physically modified easily. A naïve approach consists in plugging the resulting shape as the rest (undeformed) configuration of a shell simulator. However, once the simulation gets started, the shape irremediably sags under gravity and may greatly diverge from the desired pose, thus ruining all prior modeling efforts (see Figure 2).

Inverse cloth modeling

We believe that a new, promising way for modeling cloth both realistically and intuitively is to investigate an inverse modeling strategy. Inverse problems have been studied in a huge amount of fields ranging from Mechanical Engineering to Meteorology, Electro-Magnetism and more recently Biology. The basic principle is to retrieve unknown parameters of a model from the observation of the real object of interest. In Mechanical Engineering, inverse problems can be classified in two main categories [START_REF] Beck | Inverse problems and parameter estimation: integration of measurements and analysis[END_REF]]: (a) inverse measurement problems, where material properties (e.g., mass, stiffness, internal damping) are looked for; and (b) inverse design problems, where the undeformed configuration of the object is sought for.

For cloth, all existing studies focus on case (a), that is on the estimation of material parameters while assuming the rest shape to be perfectly flat [START_REF] Jojic | Estimating cloth draping parameters from range data[END_REF][START_REF] Bhat | Estimating cloth simulation parameters from video[END_REF][START_REF] Bouman | Estimating the material properties of fabric from video[END_REF]. Such methods, useful to identify the type of fabric involved (knit, linen, satin, etc.), are applied to simple rectangular cloth patches hung under gravity, generally in the absence of contacts or with a coarse approximation of them [START_REF] Jojic | Estimating cloth draping parameters from range data[END_REF]. Because they assume the rest pose to be flat, they cannot be used for inverting full (sewn) garments, unless the exact structure of the cloth assemblage (shape of flat panels, location of darts, etc.) can be provided [START_REF] Hasler | Reverse engineering garments[END_REF]].

In contrast, in this paper our goal is to allow the user to animate any geometrically-designed garment in a both practical and consistent way, with minimum input. That is, the user may only provide the mesh of the garment and that of the underlying body, as well as material parameters that he/she can freely choose depending on the type of motion desired. As we deliberately disregard the internal structure of the garment and want to be able to take any 3D mesh as input (not necessarily resulting from an assembly of flat patches under gravity), we consider that the rest shape of the garment cannot be determined easily and should be searched for in the space of 3D (warped) surfaces. Thus, we rather draw our attention to inverse design methods (b), where the rest shape is left as unknown.

Inverse elastic design A few approaches in Computer Graphics have recently focused on the identification of the rest pose of elastic structures such as 2D fibers [START_REF] Bertails-Descoubes | Stable inverse dynamic curves[END_REF]], mass-spring systems [Twigg and Kačić-Alesić 2011], and 3D finite element models [START_REF] Chen | An asymptotic numerical method for inverse elastic shape design[END_REF]]. Yet, these approaches are either too specific, or in contrast too generic, to take properly into account the peculiar structure of cloth models . Furthermore, none of them incorporates frictional contact in the inversion process.

The method which is closest in spirit to ours is the inverse modeling strategy recently introduced by Derouet-Jourdan and colleagues [2013] in the case of hair, which extends [START_REF] Bertails-Descoubes | Stable inverse dynamic curves[END_REF]. In this work, the authors show that it is possible to retrieve a plausible rest shape for each fiber such that at equilibrium, under gravity and frictional contact forces, the configuration of the fiber assembly matches an arbitrary set of clamped curves given as input. One key point which makes the inverse problem tractable is to rely upon curvature-based reduced models for individual fibers (e.g., super-helices [START_REF] Bertails | Super-helices for predicting the dynamics of natural hair[END_REF]), for which elastic forces depend linearly on the unknown rest configuration (natural curvature and twist). Inverse statics thus formulates as a linear equation w.r.t. the rest pose, subject to frictional constraints. Using some reasonable assumptions on the hair rest pose so as to select one solution out of the many mathematically admissible, the regularized problem boils down to a quadratic program under second-order conical constraints (SOCCQP). The latter can be solved robustly and efficiently using a solver initially devised for the forward hair dynamical problem [START_REF] Daviet | A hybrid iterative solver for robustly capturing Coulomb friction in hair dynamics[END_REF]].

To the best of our knowledge, such an inverse design strategy has never been adopted for modeling garments. Because of the nodal formulation pertaining to cloth models, strong nonlinearity is introduced in the static equations compared to [START_REF] Bertails-Descoubes | Inverse dynamic hair modeling with frictional contact[END_REF], leading to some new and important computational challenges.

Our Contribution: Inverse cloth design

We propose the first method for interpreting an arbitrary surface geometry as the mechanical result of cloth draping subject to gravity and frictional contact. Our idea is to explore inverse cloth design, where the rest position of the garment is left as an unknown of the problem, and where (unknown) frictional contact forces may also be considered to account for the input pose. Material parameters (mass, stiffness, friction coefficient) are, in contrast, not considered as unknowns and may be freely chosen by the user.

We consider a nodal cloth model similar to [START_REF] Baraff | Large steps in cloth simulation[END_REF], but with a non-flat rest shape [START_REF] Grinspun | Discrete Shells[END_REF]]. Our inverse static problem then formulates as an ill-posed nonlinear problem subject to conical constraints (Section 4). To tackle this challenging problem in a robust manner, we contribute on the following aspects:

• In the unconstrained case, we define draping as an implicit function Φ and apply the adjoint method on a modified function Φ λ so as to better drive convergence (Section 5).

• We show that frictionless bilateral constraints can be included easily to the unconstrained case (Section 6).

• In the presence of unilateral frictional contact, we show how to formulate draping as a local constrained minimization problem, and generalize the adjoint method to this constrained case (Section 7).

Finally, we successfully apply our inversion method to various cloth input geometries ranging (Section 8) and discuss limitations and further improvements (Section 9).

Cloth Inverse Static Problem

Shell cloth model To free ourselves from the knowledge of the pattern structure of garments, we choose to model a garment in one single piece. As the rest position of a sewn garment may not be flat (for instance in the case of a tight skirt or a puff sleeve), we model cloth as a thin elastic shell. This way, both the input (deformed) pose and the unknown rest pose live in the same 3D space.

Another reason to consider a non-flat rest pose is that it allows us to accommodate with possible discrepancies of the input mesh. Imagine a tablecloth, whose rest position is presumably flat. If the input pose is generated externally, there is little chance that we retrieve a perfectly flat rest pose after inversion. If we want to satisfy both the equilibrium condition and this flatness constraint, we have to tolerate a slight alteration of the input shape. However, we believe that the exact shape of the rest pose needs not be set as a hard constraint, since the rest pose is not intended to be visualized in 3D applications. We prefer to match the input pose as closely as possible, while targeting an as-flat-as-possible rest pose.

In practice we use the thin elastic shell model from [START_REF] Grinspun | Discrete Shells[END_REF]].

Notation

We represent the cloth geometry as a triangular mesh made of n nodes. Vector x t ∈ R 3n collects the 3D input (deformed) positions of the nodes, while vector x ∈ R 3n collects their unknown rest positions. If F(x, x) is a vector-valued function defined on R 3n × R 3n , then D x F(x, x) and D xF (x, x) denote its Jacobian w.r.t. the first and second variable, respectively. For a single variable function Φ( x), we shall omit the subscript in the Jacobian, that is writing DΦ( x). For a real-valued function f (x, x), we shall use the gradient notation ∇ x f (x, x) and

∇ x f (x, x) instead, with ∇ x f (x, x) = (D x f (x, x)) T ∈ R 3n .
Scalar product between two vectors x and y will be denoted by x T y.

Shell potential energy

For now, let us consider that only conservative forces act on our cloth model. Other forces such as frictional contact will be considered when formulating the statics in Section 4.2.

Conservative forces derive from the total potential energy E p (x, x) of our shell model, which can be decomposed as

E p (x, x) = E int (x, x) + E g (x)
,

where E int (x, x) is the internal elastic energy of the shell, and E g (x) its gravitational energy. Note that other external energies, such as potential fields (e.g., for modeling penalty forces), could also be considered in this formulation.

Gravitational energy Let m i be the mass of the i th node, g the constant of gravity, and e z the unitary upward vector. Gravitational energy reads

E g (x) = ∑ i∈nodes m i g x T i e z .
Note that E g is independent of the rest pose x of the cloth.

Internal elastic energy Given a deformed cloth with nodal positions x and a rest configuration x, its internal energy reads [START_REF] Grinspun | Discrete Shells[END_REF]]

E int (x, x) = E membrane (x, x) + E bending (x, x)
with the membrane (stretch and shear) energy

E membrane (x, x) = k e 2 ∑ edges (e -ē) 2 ē + k A 2 ∑ faces A -Ā 2 Ā , (1)
and the bending energy

E bending = k θ 2 ∑ edges ē h θ -θ 2 . ( 2
)
The quantities e, A, and θ represent the deformed per-edge length, per-triangle area, and per-edge dihedral angle respectively, which can all be computed from the deformed configuration x. Their rest counterparts are denoted by the barred quantities ē, Ā, and θ , which can similarly be computed from x. .

Strong nonlinearity in the rest pose

Looking at expressions (1) and (2), it is noteworthy that E int can be written formally as

E int (x, x) = 1 2 (ψ(x) -ψ( x)) T K( x)(ψ(x) -ψ( x)) = 1 2 ψ(x) -ψ( x) 2 K( x) , ( 3
)
where ψ is a vector-valued function whose components map x to either edge length e, triangle area A, or bending angle θ , and K( x) is a symmetric positive definite matrix collecting inverses of edge length or inverses of triangle area, multiplied by stiffnesses. As edge lengths, triangle areas, and bending angles are nonlinear w.r.t. positions, both ψ( x) and K( x) depend nonlinearly on x.

Because of such a nonlinearity, we anticipate here the fact that inverse statics will be a difficult problem, even in the absence of contact and friction. This will be indeed confirmed in next section.

To circumvent such a difficulty, one may think of considering a better-suited parameterization of the model. For instance, choosing ē, Ā, and θ as new unknowns of our problem, instead of x. Unfortunately, those three quantities are not independent from each other, and thus do not form a consistent parameterization .

Inverse static problem

For the sake of clarity, let us introduce the function F(x, x) = ∇ x E p (x, x). Mechanically, F represents the sum of conservative forces of the cloth.

Our goal is to interpret the input pose x t as an equilibrium configuration of the model. We would also like this equilibrium to be stable. That is, the cloth should go back to the x t configuration when slightly moved apart from x t .

Without contact In the absence of contact, x should thus satisfy

F(x t , x) = 0. (4)
Moreover, a stability condition may be written on the positiveness of D

x F(x t , x) = ∇ 2 x E p (x t , x
). Yet, such a condition is not convenient in practice because it involves the computation of eigen values of the Hessian ∇ 2

x E p (x t , x), which may be computer-expensive. Instead, we'll see in Section 5.2 how to guarantee stability of our equilibrium by formulating an appropriate draping function of the cloth. For now, we shall only consider the equilibrium condition (4).

Using Expression (3), Equation ( 4) reads

(Dψ(x t )) T K( x) (ψ(x t ) -ψ( x))) + ∇ x E g (x t ) = 0.
(5)

Equation ( 5) is, as expected, strongly nonlinear w.r.t. our unknown variable x. It is not even linear in ψ( x), because the matrix K also depends (nonlinearly) on x. Moreover, the dimension of the problem scales with the cloth size, which often reaches thousands of particles for full-size garments. Inverse statics without contact thus already turns out to be a difficult, large and nonlinear problem.

This contrasts with the simplicity of inverse static problems generated from reduced fiber models [START_REF] Bertails-Descoubes | Stable inverse dynamic curves[END_REF][START_REF] Bertails-Descoubes | Inverse dynamic hair modeling with frictional contact[END_REF]]. In such cases, K is independent of x and ψ( x) is invertible (and even equal to the identity), thus making it straightforward to retrieve x from Equation (5).

With frictional contact

In the case when the target pose x t is in contact with an external object (typically, the body of the character), we should not only rely on the rest pose x to account for the fact that x t is an equilibrium, but also on (unknown) unilateral contact forces r. Moreover, dry friction has to be considered as it clearly helps cloth remain still . Our goal is thus to find a proper balance between the role of the rest pose and that of contact and friction, so that x t closely matches a static pose of the cloth.

In the presence of frictional contact, cloth equilibrium reads

F(x t , x) = f c (x t )
where f c collects all (non-conservative) frictional contact forces, which are a priori unknown. As in [START_REF] Bertails-Descoubes | Inverse dynamic hair modeling with frictional contact[END_REF] we use the exact (implicit) Signorini-Coulomb frictional contact model introduced to graphics in [START_REF] Bertails-Descoubes | A nonsmooth Newton solver for capturing exact Coulomb friction in fiber assemblies[END_REF]. Please refer to [START_REF] Acary | Numerical methods for nonsmooth dynamical systems[END_REF] for a comprehensive analysis of this model.

Unlike conservative forces, frictional contact forces have no explicit expression, but they should satisfy some constraints. In the sequel we recall such constraints in the static case.

For the sake of simplicity, let us assume (C) that our mesh is sufficiently fine such that it is possible to consider contacts only at vertices. Then we can write f c = Hr with H a 3x3 block diagonal matrix where the i th block is defined by

H (i,i) = 0 if x i t is not a contact point R i otherwise,
with R i the rotation matrix that transforms the local (i.e., z-upward oriented) contact coordinates into world coordinates. In the static case, the cloth is necessarily sticking to the body (no relative velocity). Hence each local contact force r i must satisfy r i ∈ K µ i , where µ i is the coefficient of friction at contact i and K µ i the corresponding Coulomb friction cone .

Note that without assumption (C), H would not be block-diagonal anymore, but we would still have a linear relationship between r and f c .

To sum up, the rest position x should satisfy the following problem,

F(x t , x) = H(x t )r r i ∈ K µ i . (6)
Problem ( 6) turns out to be a nonlinear problem subject to a conical constraint. Moreover, two unknowns are now involved in the problem: the rest position x, and the local frictional contact force r. To the best of our knowledge, such a challenging (both nonlinear and nonsmooth) inverse problem has never been tackled in the past.

Inversion without Constraints

Let us start examining the unconstrained case. This case will be generalized easily in Section 6 to the case when some vertices are fixed.

Failure case: Root-finding formulation

Using the equilibrium condition (4), the inverse static problem satisfies a root-finding problem, find x such that F(x t , x) = 0. ( 7)

Relaxation of the equilibrium In the case when input data x t is slightly noisy, this problem may however not have a solution.

Instead, it is safer to consider the following minimization problem,

min x 1 2 F(x t , x) 2 g( x) . ( 8
)
Because g is smooth and coercive on R 3n , Problem ( 8) is guaranteed to possess a solution x , which is also solution to the original problem (7) if and only if the objective function g( x) vanishes at x .

We have initially considered solving this root-finding formulation, but have quickly departed from this idea for several reasons listed below.

Drawbacks of the method First, even though (7) possesses a solution x , this only guarantees that E p (x t , x ) is a critical point of energy, not necessarily a minimum of energy. Energy maxima and saddle points would thus need to be explicitly discarded at some point.

Second, solving Problem (8) (and thus Problem ( 7)) remains difficult in practice. Indeed, even in elementary 2D cloth scenarios, function g is nonconvex and features multiple local minima. Unsurprisingly, classical minimization schemes fall into such local minima, and there is no obvious way to regularize the function so as to guide minimization towards a global minimum.

Finally, in the (most probable) case when solving (8) numerically leads to some non-zero residual -which may happen either because ( 7) has no exact solution (g does not vanish at x ) or because the minimization of ( 8) gets stuck into a local minimum -the computed solution x comp turns out to be pretty useless in our case. Indeed, for such rest pose x comp , the target shape x t is not an equilibrium configuration of the cloth . The only thing we know is that this point corresponds to a low slope of energy in the x-direction.

Hence, setting x comp as the rest pose of the cloth and starting the cloth dynamic simulator from the x t state will irremediably result in some motion towards a true equilibrium configuration, which has no reason to be close to x t .

In summary

The main flaw of this root-finding formulation is to enforce the search for a vanishing point F(x, x) strictly along the x = x t subspace, even if it means finding a x that does not guarantee equilibrium at (x t , x ). Instead, we prefer to formulate the problem the other way round. That is, explore E p (x, x) on a slightly enlarged region close to but not restricted to the subspace x = x t , so as to find a pair (x , x ) which exactly satisfies equilibrium, with x as close as possible to x t .

Least-squares formulation

Instead of considering the equilibrium equation ( 7) corresponding to the balance of forces, let us model our full inverse problemincluding the stability condition -using an energy-based formulation. For a x parameter to be found, the target pose x t should match a local minimum of the function x → E p (x, x). That is, x t should correspond to a draped pose of the cloth simulator initialized with x as the rest pose.

Implicit draping function In the case when x → E p (x, x) is a strictly convex function, then a unique draped pose x exists, corresponding to the (global) minimum of the energy. We can thus define a draping function on R 3n as Φ : x → argmin x E p (x, x).

In the general case however, the potential energy is not convex. It may feature several local minima, thus precluding the definition of such a unique mapping on R 3n . This is illustrated in Figure 3. Yet, in the nonconvex case a draping function Φ can still be defined locally. Let us consider a rest pose x and a deformed pose x such that x is a strict local minimum of x → E p (x, x ). We thus have F(x , x ) = 0 and D x F(x , x ) 0. According to the implicit function theorem, there exists a unique mapping Φ from a neighborhood V of x to a neighborhood V of x such that x = Φ( x) with x ∈ V and x ∈ V . Intuitively, in our case computing Φ( x) would consist in minimizing the potential energy x → E p (x, x) in the neighborhood of (x , x ). Moreover, the theorem says that D x F(x , x ) is invertible in V × V and that Φ is differentiable in V , with an explicit expression for its Jacobian,

DΦ( x) = -(D x F(Φ( x), x)) -1 (D xF (Φ( x), x)) . (9)
Least-squares minimization Thanks to the existence of a local draping function Φ, our inverse problem literally reads find x such that x t = Φ( x), (10) and can be relaxed in the least squares sense as

min x 1 2 Φ( x) -x t 2 J( x) . ( 11
)
Similarly as before, Problem ( 11) is guaranteed to possess a solution x , which is solution to (10) if and only if the objective function J( x) vanishes at x . However, in contrast to our former problem (8), if J( x ) does not exactly vanish we still get some very useful information. Indeed, solving (11) provides us with a deformed pose Φ( x ) which is guaranteed to be a minimum configuration of the energy function x → E p (x, x) parameterized by x. Moreover, this deformed pose is as close as possible to the target x t .

Formulation ( 11) is common for inverse problems. To minimize J in practice, we use the limited-memory BFGS approach [Nocedal and Wright 2006, Section 9.1], which relies on the computation of ∇J. A classical strategy for efficiently computing ∇J is the adjoint method [Giles and Pierce 2000], which is briefly described in Section 5.3. The adjoint method has the special feature of relying only on the evaluation of Φ to compute ∇J, and not that of DΦ. However, in our case we are faced with an awkward situation, as the draping function Φ is defined locally. In Section 5.4 we shall explain how to evaluate Φ in a consistent way during subsequent steps of the minimization. Our full algorithm for minimizing J is given in Algorithm 1.

The adjoint method

To solve the minimization problem (11) robustly, it is desirable to compute the gradient ∇J accurately. Differentiating the objective function J gives

∇J = DΦ( x) T (Φ( x) -x t ) . ( 12
)
Recall that an explicit expression for DΦ( x) is provided by ( 9). However, it is dense in general and its computation requires a full matrix inversion.

Actually, we don't really need to compute DΦ( x), but only ∇J. By replacing DΦ( x) with its expression (9), the adjoint method consists in decomposing the computation of ∇J into two steps,

(D x F(Φ( x), x)) T p = x t -Φ( x) ∇J( x) = (D xF (Φ( x), x)) T p, ( 13
)
where p is called the adjoint state. Note that the adjoint method requires only one evaluation of Φ at x to compute ∇J( x).

Evaluation of the draping function

The draping problem consists in finding one local minimum of the potential energy, which amounts to solving the following problem, Given x, find x s. t. F(x, x) = 0 and D x F(x, x) 0. ( 14)

Naïve draping algorithm As already mentioned, this problem has no unique solution, and thus one cannot define a true draping function, at least globally. Still, it is easy to find an admissible x by minimizing the function x → E p (x, x) locally (using for instance the Newton method), starting from the initial guess x 0 = x t . By analogy with the draping function, let us call Φ this draping algorithm, which takes as input x and returns x.

Unfortunately, eventhough it is initialized with x t , such a method is likely to return a local minimum that is far from x t . In this case, the inverse method for solving (11) may get stuck at a wrong place and thus hardly converge to a global minimum. Furthermore, our procedure Φ is not continuous w.r.t. the x variable, as two close positions for x might lead to completely different local minima x.

Hence, during inversion it may happen that two subsequent steps of the adjoint evaluate Φ at two unrelated places. In this case the computation of ∇J becomes meaningless .

Regularized draping algorithm Our goal is to build a draping procedure that remains consistent with our inversion algorithm. That is, if x t is a local minimum of energy, the Φ procedure should be able to return it. To this aim, we penalize the energy to be minimized so as to avoid falling into a local minimum which is far from the target x t . That is, we consider the new potential energy

E λ p (x, x) = E p (x, x) + λ 2 x -x t 2 , ( 15
)
where λ ≥ 0 is a regularization factor. Choosing a high value for λ helps "convexify" the potential energy around x t . Of course it also leads to an energy that is quasi-independent of x, thus far from the original one. Conversely, setting a low value restores the original energy, but we loose the benefit of penalization. Note that if E p reaches a minimum at x t for a given x, then it is also a minimum for E λ p , ∀λ ≥ 0. This encourages use to decrease the value of λ over successive calls of Φ during inversion, when the confidence in x increases (see Algorithm 1).

Our new procedure Φ λ thus consists in minimizing x → E λ p (x, x), starting from x 0 = x t . In practice we use the Newton-CG method [Nocedal and Wright 2006, Section 6.2], to perform the minimization, that is we solve at each step k the following linear system in ∆ k+1 ,

D x F λ (x k , x) ∆ k+1 = -F λ (x k , x) (16)
where

∆ k+1 = x k+1 -x k and F λ (x, x) = ∇ x E λ p .
Compared to the pure Newton algorithm, the Newton-CG method has the advantage of always finding a descent direction ∆ k+1 , in particular in the case when D x F(x k , x) is not positive-definite. To accelerate convergence, we perform an adaptive increment along ∆ k+1 using a Wolfe linesearch. For most of our examples, we obtained convergence in a few iterations only. This good convergence rate -which contrasts with some earlier feedbacks on the Newton method [START_REF] Volino | Stop-andgo cloth draping[END_REF] -is explained by the fact that minimization starts from the target x t , which is already fairly close to the solution.

As-flat-as-possible initial guess

Multiple admissible rest poses We have seen that from one given x, multiple deformed configurations x may satisfy the draping conditions (14). Conversely, given an arbitrary shape x, multiple rest poses x may satisfy the same conditions. This non-uniqueness problem is not really an issue, as what we really want is to retrieve at least one rest pose x such that x t satisfies the draping problem ( 14). That said, when we have the choice between several admissible rest poses x, it is of course desirable to select the most plausible one.

As-flat-as-possible rest pose Because garments are designed from flat patterns, folds usually occur only because of the deformation under gravity. Thus, it is natural to target a rest pose x which is as-flat-as-possible. More precisely, since the topology of the garment is not necessarily that of a plane (in the case of a skirt or a sleeve, it has the topology of a cylinder), we would like to find the rest pose x which has the lowest possible Gaussian curvature.

Note that this choice has an influence on the realism of subsequent animation of the object, since x represent the shape that internal elasticity of the cloth wish to preserve during motion. Moreover, starting from a locally flat shape for the cloth would greatly help the texturing phase.

Flattened initial guess

In practice however, it is difficult to penalize high curvatures of x in a manner that remains consistent with our inversion algorithm. Instead, we opt for a much simpler yet effective solution: we initialize our rest pose x with a flattened version of the target pose x t , before starting our inversion algorithm.

To flatten the input configuration x t , we minimize the internal energy x → E int (x, xt ) with the rest pose xt corresponding to x t except that rest angles θ are set to zero. Thus, minimization attempts to cancel edge dihedral angles of the target x t while preserving edge lengths and face areas.

Full inversion algorithm

Let us define

J λ k ( x) = 1 2 Φ λ k ( x) -x t 2 .
Our full inversion method in the unconstrained case is summarized in Algorithm 1. To de-Algorithm 1: Robust inversion algorithm without constraints. Data:

x t , λ 0 , α ∈ ]0 ; 1[ Result: x that minimizes J // Warm start x0 ← Flatten(x t ) ; for k = 0, 1, • • • until λ k ε λ or F(Φ( xk ), xk ) ε ∇ do // Find xλ k that minimizes J λ k ( x) x0 ← xk ; for i = 0, 1, • • • until ∇J λ ( xi ) ε do // Evaluate Φ λ k at xi by draping Φ λ k ( xi ) ←Newton-CG(x → E λ k p (x, xi ), x init = x t ) ; // Evaluate ∇J λ k at xi by the adjoint ∇J λ k ( xi ) ←Adjoint(Φ λ k ( xi )) ; // Compute L-BFGS descent direction d i ←L-BFGS-Descent-Direction(∇J λ k ( xi )) ; // Update x xi+1 ← xi + d i ; // Store best minimizer of J λ k ( x) xλ k ← xi+1 ; end // Update outer iterate xk+1 ← xλ k ; // Decrease λ λ k+1 ← αλ k ;
end crease λ at each step of the outer loop, in practice we chose α = 0.5 in all our examples, except for the tablecloth for which progressive penalization was not used (one step with λ = 5 was sufficient).

Removing degrees of freedom

For some cases like a pinned piece of cloth, a flag, a tablecloth or a dress fitted at the waist, it may be useful to consider that predefined vertices of the cloth are fixed. In [START_REF] Baraff | Large steps in cloth simulation[END_REF]], any vertex i may be prevented from moving in certain space directions thanks to a simple filtering process, equivalent to an orthogonal projection of the equations [START_REF] Ascher | On the modified conjugate gradient method in cloth simulation[END_REF].

Let Π be the orthogonal projection matrix applied to the system (denoted by S in [START_REF] Ascher | On the modified conjugate gradient method in cloth simulation[END_REF]). For instance, if vertex i is fixed, we have Π i = 0, where Π i is the i th diagonal block of the matrix Π.

Modification of the inversion algorithm

To account for constrained vertices in the inversion method, it is sufficient to project draping onto the admissible space of positions, i.e., to replace (16) with

Π D x F(x k+1 -x k ) = -Π F(x k ).
As for the adjoint algorithm, it remains unchanged. . Hence, inversion with fixed vertices does not add any complexity to the unconstrained case. Since vertices in contact are assumed to be fixed due to static friction, why not use this filtering method for dealing with frictional contact ? Actually, we show in the sequel that such a method is not equivalent to the one dealing with true frictional contact, because it can yield arbitrarily oriented forces.

Bilateral constraints vs. unilateral contact

Filtering the system of equations by Π is equivalent to the adding of some frictionless bilateral constraints to the vertices. More precisely, each vertex i is subject to the linear constraint I 3 -Π i , where I 3 is the 3 × 3 identity matrix. If b i is the corresponding force (i.e., , the Lagrange multiplier) applied onto vertex i, we have

(F) i (Φ( x), x) = (I 3 -Π i ) b i ,
where (F) i (Φ( x), x) is the i th block of size 3 of F evaluated at point Φ( x), x. Thus, for a fixed vertex i, we have direct access to its applied force b i = (F) i (Φ( x), x). Figure 4a displays these bilateral forces (black arrows) for a simple 2D case. In this example, our target equilibrium shape x t is depicted in red. It corresponds to a curve, slightly curly at both ends, lying on a circular object. To invert this pose, we started from the flattened shape represented as the blue dotted line, and fixed all contacting point during inversion. Inversion retrieves the solid blue configuration x , straight in the middle and curly at both ends. Note that the resulting equilibrium pose Φ( x ) (in green) is identical to the target x t .

At a first glance, this rest pose looks consistent with the input data. Indeed, only points outside the contact zone have been displaced from the initial flat guess. However, it is noteworthy that the "gluing" forces, represented with black arrows, cannot be interpreted as consistent contact forces. Some of them are penetrating the circular obstacle while others, correctly pointing outwards the obstacle, share the same orientation upon the circle as if the friction coefficient was varying along the surface.

In contrast, our inversion with frictional contact, displayed in Figure 4b, automatically yields forces belonging to the local Coulomb friction cone (see next section). Because this method does not tolerate forces to have any arbitrary orientation like in the previous case, it retrieves a more plausible rest pose x . In particular, it was able to reverse the stretching applying to the contacting zone of the target, and to retrieve a slightly more compressed rest pose than the filtering method.

The filtering method is thus limited to constraints that should still hold in the rest configuration. Yet, it can still be used to generate a good warm start for inverting the system with contact and friction (see next section).

Including Frictional Contact

Starting from the equilibrium condition (6), our inverse static problem formulates as find ( x, r) such that

F(x t , x) = H(x t )r r i ∈ K µ i . ( 17
)
Obviously this problem has no unique solution, but as previously, we show that an initial as-flat-as-possible warm start greatly helps select a plausible rest pose x.

Similarly as in Section 5, we do not want to solve strictly Problem ( 17), but rather relax the problem using a least-squares formulation. This raises two major challenges:

How to compute draping ? It is noteworthy that the balance of forces under frictional contact (17) involves two dual variables: (a) the primal (position) x variable, and (b) the dual (force) r variable.

As is, such a mixed formulation is not adapted for computing a static pose x given a rest pose x. In Section 7.1 we show how to formulate draping as a constrained minimization problem on the primal variable x only, using the fact that contact point locations are known from the target pose x t .

How to generalize the adjoint method ? Recall that the adjoint method is defined starting from an implicit equation F(x, x) = 0 where F is assumed to be differentiable. When unilateral contact with friction is considered, we come up with a new implicit equation F c (x, x) accompanied with a conical constraint (see Problem (17). This new format thus does not fit in with the adjoint framework straightforwardly. In Section 7.3 we introduce a nonsmooth projection to convert our draping problem into a zero-finding formulation, and show that the adjoint keeps on working well in this nonsmooth case.

Towards a least-squares formulation

Primal vs. dual formulation The static formulation ( 17) is pretty inconvenient, as the constraint is expressed on the dual variable r. Moreover, E p is nonlinear with respect to x, making it difficult to express x as a function of r. Note that this contrasts with reduced fiber models that feature linear elastic forces, which allows for a straightforward elimination of the primal variable x in the system [START_REF] Bertails-Descoubes | Inverse dynamic hair modeling with frictional contact[END_REF].

One solution would be to linearize F locally, so as to formulate a local problem depending only on the dual variable r, as done for instance in [START_REF] Kaufman | Adaptive nonlinearity for collisions in complex rod assemblies[END_REF]. However, such a method would be computationally expensive in our case, as it would yield a dense Delassus operator.

Instead, by using some tools of convex analysis, we demonstrate here that a much simpler formulation of (17) can be derived, which only depends on the primal variable x.

A new constrained minimization problem Let us assume for now that x → E p (x, x) is a convex function. Let C be a convex set of R 3n . From convex analysis [Hiriart-Urruty and Lemaréchal 1993, Theorem 2.1.4] we have the following result,

min x∈C E p ⇐⇒ -F(x) ∈ N C (x), (18)
where N C (x) is the normal cone of C at x, defined as

N C (x) = {d ∈ R 3n , d T (c -x) ≤ 0 ∀c ∈ C}.
In the particular case when the convex set C is a cone K 1 µ of aperture 1 µ , the normal cone boils down to

N K 1 µ (x) = -K µ ∩ x ⊥ where x ⊥ = {y ∈ R 3n , y T x = 0}.
Thus, the equivalence (18) can be expressed as min

x∈K 1 µ E p (x) ⇐⇒ K µ F(x) ⊥ x ∈ K 1 µ . ( 19
)
This condition is illustrated in .

Let us assume our target pose x t has one contact located at the i th node x t,i . From (17) we get

(F(x t , x)) i ∈ R i K µ i , (20)
where R i K µ i is still a cone of aperture µ i , corresponding to a rotated cone K µ i .

Besides, let us consider the following minimization problem

min x∈C E p (x, x) with C = x ∈ R 3n , H T (x -x t ) ∈ K 1 µ , (21)
where

K 1 µ = K 1 µ 1 × ... × K 1 µ n . We have x ∈ C ⇐⇒ H T (x -x t ) ∈ K 1 µ ⇐⇒ R T i (x i -x t,i ) ∈ K 1 µ i ∀i in contact ⇐⇒ (x i -x t,i ) ∈ R i K 1 µ i ∀i in contact.
Thus, using Equivalence (19), it follows that Problem ( 21) is equivalent to solving

∀i in contact, R i K µ i (F(x, x)) i ⊥ R i (x i -x t,i ) ∈ K 1 µ i ∀i not in contact, (F(x, x)) i = 0.
For contact-free points, the initial problem ( 17) is automatically satisfied when x t = x. For points in contact, when x = x t , the right-hand side of the above complementarity condition vanishes, leading exactly to (20) and thus, again, to our initial problem (17

Evaluation of the draping function

Due to the simple nature of the constraints appearing in ( 21), the draping problem can be solved with a simple algorithm such as a (conjugate) projected gradient. Indeed, it suffices to remark that the orthogonal projection of x on the constraints is explicitely given by

Π C x i = x i if x i t is not a contact point x i t + Π K 1 µ i R i (x i -x i t ) otherwise, where Π K 1 µ i (•) is the orthogonal projection on cone K 1 µ .
The draping procedure is provided in Algorithm 2.

Algorithm 2: Projected gradient for evaluating the draping with contact and friction.

Data: x, H, K 1 µ , α ∈ [0, 1] Result: x that satisfies Equation (17) x 0 ← x t ; for k = 0, 1, • • • do d k ← F(x k , x) ; if x k -Π C (x k -d k ) sufficiently small then return x k ; end ρ ← 1; repeat x k+1 ← Π C (x k -ρd k ); ρ ← αρ ; until E p (x k+1 , x t ) < E p (x k , x t ); end

Fitting in with the adjoint method

In order to use the formalism of the adjoint method, we can first rewrite Equation (17) using the orthogonal projection operator

Π H,K µ x i := R i Π K 1 µ i R i x i .
Equation ( 17) is then equivalent to the root-finding problem

h(x, x) := (I -Π H,K µ (•)) • F(x, x) = 0. ( 22
)
The adjoint algorithm can be now written in a similar fashion than without contacts. That is, for a given value of x, 1. Compute x = Φ( x) ;

2. Compute the adjoint state p as ∂ h ∂ x (x , x) p = x tx ;

3. Compute the gradient of J as ∇J( x) = ∂ h ∂ x (x , x) p. Similarly as in the bilateral case, it can be shown that the adjoint method does not need to be modified. . Furthermore, in practice we have not observed any issue related to the fact that h is locally nonsmooth.

Results

Framework

We have tested our inversion method on three main 3D geometries:

• A tablecloth, generated from simulation. Contacts (onto a centered round table) are modeled by fixing positions, using the method of Section 6;

• A puff sleeve, modeled by an artist. Contacts (localized at both extremities of the sleeve) are also modeled by fixing positions;

• A skirt, modeled by an artist. Unilateral contacts with dry friction are considered between the garment and the body, using the method of Section 7.

Tablecloth:

Skirt: 5208 vertices and 755 contact points Puff sleeve: 2045 vertices and 102 contact points

For each of these geometries, we evaluated the correctness of our method in two ways:

• First numerically, by evaluating the closeness of the equilibrium to the target pose. These results are compiled on Figure 5.

• Then visually: we put each equilibirum in a forward cloth simulator to ensure that (a) when the simulation starts, the cloth does not sags, and (b) when we add a new external force, such as wind, the cloth behaves in plausible manner. These simulations are shown in the accompanying video.

For the puff sleeve and the skirt, we were not able to match the target perfectly. However, in those cases we were still able to match the intended look of the cloth, and found a much better starting point for the simulation than originally available.

Performance While devising the fastest algorithm was not our primary goal, we are aware that being too slow could significantly limit the applicability of our method. The geometry which took the most time to invert was the skirt, which features about 5200 vertices, with a wall-clock time of just under three hours.

Discussion

Modeling cloth as a thin elastic shell In the case when cloth is modeled manually using a 3D free form interface, there is no guarantee that the result satisfies realistic cloth constraints. Typically, Advantages In comparison with the traditional work flow of pattern design, where the user continuously transitions between pattern editing and draping simulation, our modeling process is very different, as it is split into two successive phases: (a) the geometric modeling phase, which is fully interactive or automatic, and (b) our automatic (and fast) conversion into a physics-based cloth model. Unlike physics-based tailoring, our approach is geometry centered and physics is transparent to the user: geometry is freely explored by the user, while physics should obey (and not the reverse). Any shape is feasible, but of course physical realism may be altered (cloth may be made with shell patterns).

Limitations For the sake of simplicity, we assumed here the garment material properties (mass, stiffness) to be known. This simplifying assumption appears reasonable to us in a first step. Indeed, whereas a designer cannot easily guess the rest pose of a deformed configuration, she/he generally has some idea about the kind of material (silk, linen) she/he would like to simulate. Moreover, as suggested by some studies in fiber design [START_REF] Bertails-Descoubes | Stable inverse dynamic curves[END_REF]] and as we show here, for a well-chosen rest shape and regardless of the material properties, any given configuration may be interpreted as a static pose. This means that even if the designer chooses a "wrong" material, our method will still provide a rest pose which allows to interpret the input as a static pose. But of course, in such a case the rest pose may be pretty unrealistic . .

In some cases (like the puff sleeve) where the target is mainly explained by bending forces, our algorithm has more difficulty to converge to an equilibrium that is very close to the target. In practice we have noticed that the success of our convergence is highly dependent upon the formulation of the bending energy. We would like to investigate this issue in the future.

Still, our method produces results which greatly improve naïve ap-proaches consisting in (a) initializing x with x t , or (b) initializing x with a flattened version of x t (see Figure 5).

Future work So far we have not considered self-contacting for inversion. Adding self-contacts would not bring any difficulty to our theoretical framework. However, we would have to investigate a faster draping function so as to be able to scale up the number of contacts (currently limited to a few thousands).

Another interesting venue of our approach would be to generate realistic patterns from our retrieved rest pose. Indeed, current methods compute such patterns by flattening the deformed shape (the target), thus neglecting deformations due to stretching and shearing [START_REF] Brouet | Design preserving garment transfer[END_REF]. With our approach we could generate more realistic patterns, and thus fully resolve the 3D→2D computation. In that sense our approach could be served to enrich current pattern design methods such as [START_REF] Umetani | Sensitive couture for interactive garment editing and modeling[END_REF]

Figure 2 :

 2 Figure 2: Naïve initialization of the rest pose by the target (left) yields sagging when statics is simulated (right), thus diverging from the target pose. In contrast, our inverse cloth design strategy allows to recover a plausible rest pose, under which the target pose matches a static configuration of the cloth simulator.

( a )

 a Rest pose x (b) Two draped poses, side view (c) Two draped poses, bottom view

Figure 3 :

 3 Figure 3: For this tablecloth, at least two different local minima of potential energy (one in blue and one in orange) may be obtained from the same rest pose x (in light grey).

  Inversion with contact modeled by fixing positions. Equilibrium (target) Rest (inversion) Rest (initial value) Equilibrium (inversion) (b) Inversion with unilateral contact and dry friction.

Figure 4 :

 4 Figure 4: Comparison of the yielded contact forces in the bilateral case and in the unilateral case with dry friction.

Figure 5 :

 5 Figure 5: Evaluation of our inversion method on our 3 examples. Left column shows the target pose; middle column compares the sagged (naïve) pose against the target, and right column compares our generated equilibrium Φ( x ) against the target. Blue color indicates a very low error , whereas red indicates an error of the order of 10% of a typical scene length.

  ). Very interestingly, this means that ( x, r) is solution to (17) if and only if x t is solution to (21). In the general (nonconvex) case, we still benefit from the following sufficient condition: if x t is solution to (21), then ( x, r) is solution to (17). We thus transform the solving of (17) into the searching of x such that x t matches a local minimum of Problem (21). Locally solving (21) thus defines a new "draping" procedure Φ( x) which takes into account frictional contact. Similarly as in the unconstrained case, we can thus imagine to derive a least-squares formulation for finding a x which allows x t to match Φ( x). Still, two questions remain: (a) How to evaluate Φ( x) in practice ? and (b) How to modify the adjoint method so that our new draping function can fit in ? These two issues are examined in the following sections.

	Draping with unilateral contact and friction

Acknowledgements

The authors would like to thank Laurence Boissieux for the design of the three cloth meshes used as targets (tablecloth, puff sleeve, skirt), as well as for the final rendering of the different scenes.


Pierre Alain
email: pierre.alain@univ-rennes1.fr

Jonathan Chevelu
email: jonathan.chevelu@irisa.fr

David Guennec
email: david.guennec@irisa.fr

Gwénolé Lecorvé
email: gwenole.lecorve@irisa.fr

Damien Lolive
email: damien.lolive@irisa.fr

The IRISA Text-To-Speech System for the Blizzard Challenge 2016

Keywords:

This paper describes the implementation of the IRISA unit selection-based TTS system for our participation in the Blizzard Challenge 2016. We describe the process followed to build the voices from given data and the architecture of our system. The search is based on a A* algorithm with preselection filters used to reduce the search space. A penalty is introduced in the concatenation cost to block some concatenations based on their phonological class. Moreover, a fuzzy function is used to relax this penalty based on the concatenation quality with respect to the cost distribution.

Introduction

In recent years, research in text-to-speech synthesis essentially focused on two major approaches. The first one is the parametric approach, for which HTS [START_REF] Yamagishi | Robustness of HMM-based speech synthesis[END_REF] and DNN-based systems [START_REF] Ze | Statistical parametric speech synthesis using deep neural networks[END_REF] are now dominating the academic research in recent years. This method offers advanced control on the signal and produces very intelligible speech but with a low naturalness. The second approach, unit selection, is a refinement of concatenative synthesis [START_REF] Sagisaka | Speech synthesis by rule using an optimal selection of non-uniform synthesis units[END_REF][START_REF] Black | CHATR: a generic speech synthesis system[END_REF][START_REF] Hunt | Unit selection in a concatenative speech synthesis system using a large speech database[END_REF][START_REF] Taylor | The architecture of the Festival speech synthesis system[END_REF][START_REF] Breen | Non-uniform unit selection and the similarity metric within BTs Laureate TTS system[END_REF][START_REF] Clark | Multisyn: Open-domain unit selection for the Festival speech synthesis system[END_REF][START_REF] Patil | A syllable-based framework for unit selection synthesis in 13 indian languages[END_REF]. Speech synthesized with this method features high naturalness and its prosodic quality is unmatched by other methods, as it basically concatenates speech actually produced by a human being.

The challenge for 2016 is to build an expressive voice using children's audiobooks in English. The main difficulty with audiobooks, and in particular for children, is the change of characters and here the imitation of animals (i.e. roars) as well as other sounds (i.e. bell ringings) that may occur. For instance, in sample data provided, a signal is given to tell the child that he/she has to turn the page. Considering the expressivity of the voice and the different sounds and characters we can find in such books, the main challenge is phone segmentation and expressivity control.

In this paper we present the unit-selection based IRISA system for the Blizzard Challenge 2016. Basically, the system is based on preselection filters to reduce the acoustic unit space to explore and on an A* algorithm to find the best unit sequence. The objective function minimized by the algorithm is composed of a target cost and a join cost. The join cost relies mainly on acoustic features to evaluate the level of spectral resemblance between two voice stimuli, on and around the position of concatenation. For instance, distances based on MFCC coefficients and especially F0 are used [START_REF] Stylianou | Perceptual and objective detection of discontinuities in concatenative speech synthesis[END_REF][START_REF] Tihelka | Modelling F0 Dynamics in Unit Selection Based Speech Synthesis[END_REF]. In particular, for the challenge, we have introduced a penalty on units whose concatenation is considered as risky. This follows the work of [START_REF] Yi | Natural-sounding speech synthesis using variable-length units[END_REF][START_REF] Cadic | Vocalic sandwich, a unit designed for unit selection TTS[END_REF] which showed that artefacts occur more often on some phonemes than others. For this purpose, we define a set of phoneme classes according to their "resistance" to concatenation. A phoneme is called resistent if the phones of its class are usually unlikely to produce artefacts when concatenated. This approach has been originally proposed in the context of recording script construction in [START_REF] Cadic | Vocalic sandwich, a unit designed for unit selection TTS[END_REF] to favor the covering of what has been called "vocalic sandwiches".

The remainder of the paper is organized as follows. Section 2 describes the voice creation process from the given data. Section 3 details the TTS system. Section 4 presents the evaluation and results.

General voice creation process

This year, the challenge focuses on audiobook reading for children in English. The goal is then to build a voice based on the 5 hours of speech data provided as a set of wave files with the corresponding text. The very first step has been to clean the text and make sure that it was corresponding to the speech uttered.

Then to build the voice, we first phonetized the text thanks to a grapheme-to-phoneme converter (G2P) and then, using another tool, automatically segmented speech signals according to the resulting expected phonemes. As for the G2P tool, we used eSpeak [START_REF] Duddington | espeak text to speech[END_REF]. Once phonetized, speech signals have been segmented using the language independent segmenter MAUS [START_REF] Schiel | Automatic phonetic transcription of non-prompted speech[END_REF]. We have also used the ROOTS toolkit [START_REF] Chevelu | ROOTS: a toolkit for easy, fast and consistent processing of large sequential annotated data collections[END_REF] to store all the necessary information and produce conversions from IPA (output from eSpeak) to the SAMPA alphabet (used by MAUS).

3. The IRISA system 3.1. General architecture The IRISA TTS system [START_REF] Guennec | Unit Selection Cost Function Exploration Using an A* based Text-to-Speech System[END_REF], used for the experiments presented in this paper, relies on a unit selection approach with an optimal graph-search algorithm (here an A* algorithm). The optimization function is divided, as usually done, in two distinct parts; a target and a concatenation cost [START_REF] Black | CHATR: a generic speech synthesis system[END_REF] as described below:

U * = argmin U (Wtc card(U ) n=1 wnCt(un) + Wcc card(U ) n=2 vnCc(un-1, un)) (1)
where U * is the best unit sequence according to the cost function and un the candidate unit trying to match the n th target unit in the candidate sequence U . Ct(un) is the target cost and Cc(un-1, un) is the concatenation cost. Wtc, Wcc, wn and vn are weights for adjusting magnitude for the parameters. Subcosts are weighted in order to compensate magnitudes of all sub-costs as in [START_REF] Blouin | Concatenation cost calculation and optimisation for unit selection in TTS[END_REF]. In practice, the weight for each sub-cost c is set to 1/µc, where µc is the mean sub-cost c for all units in the TTS corpus. The problem of tuning these weights is complex and no consensus on the method has emerged yet. [START_REF] Alías | Efficient and reliable perceptual weight tuning for unit-selection text-to-speech synthesis based on active interactive genetic algorithms: A proof-ofconcept[END_REF] is a good review of the most common methods.

Join cost

The concatenation cost Cc(u, v) between units u and v is composed of MFCCs (excluding ∆ and ∆∆ coefficients), amplitude, F0, syllable speech rate and syllable F0 level euclidean distances, as below:

Cc(u, v) = C mf cc (u, v) + Camp(u, v) + CF 0(u, v) +Crate(u, v) + C lev (u, v) + C dia (u, v),
where C mf cc (u, v), Camp(u, v), CF 0(u, v), Crate(u, v), C lev (u, v) and C dia (u, v) are the sub-costs, resp., for MFCC, amplitude, F0, speech rate, F0 level and dialog section. The speech rate and the F0 mean level are computed on a syllable basis and are averaged on a window of ±1 syllable around the current phoneme. The dialog sub-cost is a penalty that is added if the phonemes u and v are taken from inconsistent parts of the corpus with respect to the narrative style.

Target cost

For candidate units, we compute a numerical target cost as a weighted sum of the features given in table 1. The features used in the computation are nearly the same as the ones used for preselection as explained in section 3.4. The weights Wtc and Wcc used in (1) are arbitrarily set to balance the importance of the join cost compared to the target cost.

Preselection filters

When exploring new units in the graph, the algorithm accesses to the corpus via an ordered list of preselection filters, where the role of each filter is to reject speech units which do not respect a given specific property. Their purpose is twofold. First, it considerably prunes the graph explored by the unit selection algorithm, making the selection process faster. Second, it serves as a set of binary target cost functions relying on the assumption that if a unit does not respect the required set of features, it can't be used for selection. The preselection filters should therefore be seen as part of the cost for a unit. In our system, when no corpus unit (or not enough units) respects a given set of preselection filters, the set is temporarily relaxed (removing one by one the features that seem the less helpful) until units are found. This mechanism ensures finding a path in all cases under the assumption that the whole corpus contains at least one instance of the most basic units, i.e. diphonemes. In our case, the threshold number of units is set to 100.

In case a diphoneme is not present in the corpus, a fallback mechanism has been implemented. Precisely, the requested diphoneme is built artificially by concatenating two phonemes in context of a pause. As it does not take into account coarticulation effects, the result is not excellent but at least, it enables to produce speech. The set of preselection filters we use in this work is the following:

1. Unit label (mandatory). 2. Is the unit a pause (mandatory)? 3. Is the phone nasal ? 4. Is the phone long ? 5. Is the phone stressed (primary stress) ? 6. Is the phone stressed (secondary stress) ? 7. Is the phone in a dialog part ? 8. Is the phone in the last syllable of its sentence? 9. Is the phone in the last syllable of its breath group? 10. Is the phone in a syllable with a rising intonation ? 11. Is the phone in a syllable with a falling intonation ? 12. Is the current syllable in word end?

The two first filters, written as mandatory, cannot be relaxed as they represent the minimal information to retrieve units.

Fuzzy concatenation cost

Analysis of synthesized sentences containing artefacts shows that concatenation on some phonemes, especially vowels and semi-vowels, is more likely to engender artefacts than others (plosives and fricatives for example, especially unvoiced ones) [START_REF] Yi | Natural-sounding speech synthesis using variable-length units[END_REF]. Phonemes featuring voicing, high acoustic energy or important context dependency are generally subject to more distortions. Based on this ascertainment, [START_REF] Cadic | Vocalic sandwich, a unit designed for unit selection TTS[END_REF][START_REF] Cadic | High Quality TTS Voices Within One Day[END_REF] proposed a corpus covering criterion where the objective is to get a maximum covering of "sandwich units". A sandwich unit is a sequence of phonemes where one or several syllabic nuclei are surrounded by two phonemes considered as robust to concatenation artefacts. Concerning unit selection concatenation costs, a few work can also be cited, for example [START_REF] Donovan | A new distance measure for costing spectral discontinuities in concatenative speech synthesizers[END_REF][START_REF] Yi | Natural-sounding speech synthesis using variable-length units[END_REF], but in these works, costs and penalties are not flexible enough. In unit selection, too many constraints generally means loss of quality (e.g. too many preselection filters is to prevent).

In the approach we introduced in [START_REF] Alain | The IRISA Text-To-Speech System for the Blizzard Challenge 2015[END_REF][START_REF] Guennec | On the suitability of vocalic sandwiches in a corpus-based tts engine[END_REF], we have defined a fuzzy concatenation cost taking into account three phonetic clusters:

V (vowel) : Vowels, on which concatenation is hardly acceptable.

A (acceptable) : Semi-vowels, liquids, nasals, voiced fricatives and schwa. These units are viewed as acceptable concatenation points, but still precarious.

R (resistant) :

The remaining phonemes (unvoiced consonants, voiced plosives), where concatenation is definitely possible.

First, we give a fixed penalty to each phoneme class: 0 for phonemes in R, a penalty slightly higher than the highest value Cc observed in the corpus for all phonemes in A. Vowels (V) are given a huge penalty, big enough to prevent compensation by other costs in the candidate sequence. It corresponds to a penalization of candidate units based on the phonemes on which concatenation may be performed if choosing this unit. In this case, a new concatenation cost function C c is formulated as:

C c (u, v) = Cc(u, v) + K(u, v), (2)
where K(u, v) = p(v) is the penalty depending on the phoneme that begins the unit v as described before.

In order to relax this penalty when a concatenation between two candidate units is statistically among the best ones, we introduce a fuzzy weighting function, ranging from 0 to 1. It describes how much the unit belongs to one of the clusters defined earlier.

Assuming MFCC, Amplitude and F0 cost distributions follow normal distributions, we define two thresholds for each subcost. For instance, the two thresholds T 1 F 0 and T 2 F 0 for the F0 sub-cost may be defined as:

T 1 F 0 = µC F 0 -σC F 0 (3)
T 2 F 0 = µC F 0 + σC F 0 (4)
Formally, the fuzzy function is defined, for the F0 sub-cost:

fF 0 (u, v) =        0 if CF 0 (u, v) < T 1 F 0 , 1 if CF 0 (u, v) > T 2 F 0 , 1.0 - (T 2 F 0 -C F 0 (u,v)) (T 2 F 0 -T 1 F 0 ) otherwise.
(5) The choice for that tolerance interval is motivated by the observation of real cost distributions. To be complete, the choice of the thresholds should be differentiated depending on the type of sub-cost and optimized separately. In this paper, we used thresholds corresponding to 15% of the distribution for each sub-cost.

The penalty is then modified in the following way:

K(u, v) = (f mf cc (u, v) + famp(u, v) + fF 0(u, v)) * p(v)
where f mf cc (u, v), famp(u, v) and fF 0(u, v) correspond to the fuzzy function of the form described previsouly respectively for MFCC, amplitude and F0. The value p(v) is still the generic penalty value that depends on the phoneme class and which is not weighted. With this fuzzy function, the main idea is to decrease the penalty when the unit has a concatenation sub-cost value which is statistically among the best ones. The sub-cost distributions are estimated from the voice corpus by computing concatenation sub-costs for F0, amplitude and MFCC using all the units in the corpus.

Break prediction

As we had no pause prediction module until now and as it is very important for audiobooks, specifically for children, we tried to introduce a simple prediction block in our system. We have chosen to formulate it as a classification task by assigning a label to each word telling if it is followed by a pause or not, and also the type of the pause. To simplify, we make 4 clusters:

• No Break (NB): length less than 120ms; • Short Break (SB): length between 120ms and 250ms;

• Medium Break (MB): length between 250ms and 550ms;

• Long Break (LB): length greater than 550ms.

To predict those labels, we used a random forest classifier [START_REF] Breiman | Random forests[END_REF] with 100 estimators and 77 features as input which consisted of linguistic features, positional features and also a narrative vs. dialog feature with a window of ±2. We used the Scikit toolkit [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] to learn the classifier.

Using the data provided, the overall accuracy score is 84%. The confusion matrix is given in table 2. As we can observe, NB class is very well predicted with nearly no confusion which is understandable as it is the dominating value in the data set. The most confused class is the SB one. Interestingly, the longer the break, the less confusion. Despite the unbalanced classes, the random forest is doing quite well.

Evaluation and results

The evaluation assessed a number of criteria (overall impression, speech pauses, intonation, stress, emotion, pleasantness and listening effort) for book paragraphs as well as similarity to the original speaker, naturalness and intelligibility. The evaluation has been conducted for different groups of listeners: paid listeners, speech experts, and volunteers. In this section, we only give results including all participants. In every figure, results for all 17 systems are given. Among the system, we have the following : system A is natural speech, system B is the Festival benchmark system (standard unit selection), system C is the HTS benchmark and system D is a DNN benchmark. System H (in orange on all figures) is the system presented by IRISA.

Evaluation with paragraphs

Results are shown on figures from 1 to 7. For each criterion, our system achieves average results, except for the speech pauses evaluation (fig. 3) showing a low score. After the submission of speech samples for the evaluation, we have detected a bug in the break prediction module causing an inconsistency between results during the learning step and its use in the synthesis chain. Concerning the other criteria, these average results are likely to be explained by the quality of the phone segmentation. A manual checking revealed an alignment problem that couldn't be solved before the submission.

Similarity to original speaker and naturalness

The similarity of the speech produced, as shown on figure 8, is among the best ones with a mean score of 3.6 and the median value at 4. Similarly, naturalness is also quite good as shown on figure 9 with an average of 2.8 and a median of 3. Those results are coherent with the nature of our system. Unit selection uses natural stimuli and consequently preserves the original speaker's voice.

Intelligibility

Concerning intelligibility, our system performed poorly compared to the majority of other systems with an average word error rate of 52%. Detailed results are given on figure10. The main explanation is the low quality of phone segmentation we obtained.

Conclusion

We described the unit-selection based IRISA system for the Blizzard challenge 2016. The unit selection method is based on a classic concatenation cost to which we add a fuzzy penalty that depends on phonological features. In order to improve the system, we added specific costs to deal with speech rate, melody and speaking style (narrative and dialog) consistency. Despite the improvements we've made, our system obtained average results. The main explanation is that the phone segmentation system we used performed badly on the given data. This was the cause of a drop in nearly all criteria. Concerning pause prediction, we have found a bug in the pause prediction module that will be fixed for next year edition.

Figure 1 :

 1 Figure 1: Mean Opinion Scores, overall impression evaluation, all listeners (from "bad" to "excellent").

Figure 2 :

 2 Figure 2: Mean Opinion Scores, pleasantness evaluation, all listeners (from "very unpleasant" to "very pleasant").

Figure 3 :

 3 Figure3: Mean Opinion Scores, speech pauses evaluation, all listeners (from "speech pauses confusing/unpleasant" to "speech pauses appropriate/pleasant").

Figure 4 :

 4 Figure4: Mean Opinion Scores, intonation evaluation, all listeners (from "melody did not fit the sentence type" to "melody fitted the sentence type").

Figure 5 :

 5 Figure 5: Mean Opinion Scores, stress evaluation, all listeners (from "stress unnatural/confusing" to "stress natural").

Figure 6 :

 6 Figure6: Mean Opinion Scores, emotion evaluation, all listeners (from "no expression of emotions" to "authentic expression of emotions").

Figure 7 :

 7 Figure 7: Mean Opinion Scores, listening effort evaluation, all listeners (from "very exhausting" to "very easy").

Figure 8 :

 8 Figure 8: Mean Opinion Scores, similarity to the original speaker evaluation, all listeners.

Figure 9 :

 9 Figure 9: Mean Opinion Scores, naturalness evaluation, all listeners.

Figure 10 :

 10 Figure 10: Word Error Rates, naturalness evaluation, all listeners.

Table 1 :

 1 List of features used in the target cost

	Text related features:
	TEXT DIALOG
	Phoneme position:
	LAST OF BREATHGROUP
	LAST OF WORD	LAST OF SENTENCE
	IN CODA	IN ONSET
	SYLLABLE BEGIN	SYLLABLE END
	WORD BEGIN	WORD END
	Phonological features:
	LONG	NASAL
	LOW STRESS	HIGH STRESS
	Syllable related features:
	HAS CODA
	LAST SYL OF SENTENCE
	LAST SYL OF BREATHGROUP
	SYLLABLE RISING SYLLABLE FALLING

Table 2 :

 2 Normalized confusion matrix

		NB	SB	MD	LB
	NB	98.1	0.6	1.0	0.2
	SB	28.8 28.4 32.9	9.9
	MB 14.7 19.2 46.1 20.1
	LB	6.9	6.0	31.2 55.9

Acknowledgements

This study has been partially funded thanks to the ANR (French National Research Agency) project SynPaFlex ANR-15-CE23-0015.