Marion Gstalter

Grégoire Chabrol

Armel Bahouka

Kokou-Dodzi Dorkenoo

Jean-Luc Rehspringer

Sylvain And

Lecler

Long focal length high repetition rate femtosecond laser glass welding

A long focal length focusing device is proposed for the process of glass welding by femtosecond laser pulses at high repetition rate and report on the significant advantages. The study is performed using a 100 mm focusing length F-theta lens. The results are compared to those obtained with high numerical aperture microscope objective. The long focal length with the associated Rayleigh length method allows a robust high process speed: welding at 1000 mm/s has been achieved, several order of magnitude larger compared to what was reported till now. Moreover, the heat accumulation process on a larger laser spot leads to a lower temperature increase after each pulse and thus a lower thermic gradient. As a result, the residual stress in the welding seams is reduced, preventing the formation of fractures inside the seams: mechanical resistance at 30 MPa has been measured

Introduction

Glass welding by femtosecond laser pulses is a recent method which has been developed to overcome the weaknesses of the conventional glass bonding technics [START_REF] Richter | Laser welding of glasses at high repetition rates -Fundamentals and prospects[END_REF][START_REF] Okamoto | Evaluation of molten zone in micro-welding of glass by picosecond pulsed laser[END_REF]. Classical processes such as adhesive bonding or fusing require either the addition of an adhesive material or high temperature. The first one limits the chemical and thermal resistance of the bonding and the second one is unsuitable for temperature sensitive material [START_REF] Hülsenberg | Microstructuring of glasses[END_REF]. By comparison, glass welding by ultrashort laser pulses presents many advantages in term of mechanical, thermal and chemical resistance, control and precision, biocompatibility and process speed [START_REF] Zimmermann | Ultrastable bonding of glass with femtosecond laser bursts[END_REF][START_REF] Carter | Picosecond laser micro-welded similar and dissimilar material[END_REF][START_REF] Cvecek | Gap bridging in joining of glass using ultra short laser pulses[END_REF]. This method is particularly suitable for micro-welding applications, such as micro-optics, microfluidics or MEMS packaging [START_REF] Tamaki | Welding of transparent materials using femtosecond laser pulses[END_REF].

Femtosecond laser glass welding relies on the irradiation of a volume at the interface of two glass plates in close contact with a focused laser beam, as illustrated in Fig. 1. Due to the low linear absorption at 1030 nm in our case, the beam can propagate through the first glass plate until high optical intensity is reached in the focused spot, generating nonlinear absorption of the laser energy [START_REF] Tamaki | Welding of transparent materials using femtosecond laser pulses[END_REF]. The use of a high repetition rate laser (above 300 kHz for borosilicate glass) introduces thermal accumulation effects leading to a localized temperature increase of the glass in the focusing point up to its melting point [START_REF] Richter | Laser welding of glasses at high repetition rates -Fundamentals and prospects[END_REF]. The bonding is obtained during the fast cooling of the melting pool.

Probably due to the former difficulty in reaching the non-linear absorption threshold with "classic" lasers (pulse duration ≥ ps), femtosecond laser glass bonding at high repetition rate has mainly been demonstrated using microscope objective with high numerical aperture (Fig. 2) [START_REF] Zhang | Femtosecond laser Bessel beam welding of transparent to nontransparent materials with large focal-position tolerant zone[END_REF]. The work presented in this paper demonstrates the benefits of using a scanner head containing a long focusing length lens with low numerical aperture (Fig. 3). At first glance, the welding process can be described similarly in both configurations: nonlinear absorption of the laser beam and thermal accumulation effect due to high repetition rate. The

Experim

Please see the

Material

Experiments h

The glass pla and a flatnes without the u µm between fringes. glass welding pro microscope objectiv of using a long high numeri ng quality (lo will be presen ummarizes all µm thin boros ity with an ari ality is suitabl device [START_REF] Cvecek | Sample preparation method for glass welding by ultrashort laser pulses yields higher seam strength[END_REF]. Ho e locally obser much higher (a ad.

uency involved roscope objecti

Laser sy

The experime developed spe high repetitio generating 30 adapted using The focusing mirrors, and c diameter of th of 2 µm and c

Results

High wel

As can be see brings one ma scanning spee interface. Na available, alm seams compa and energy, a measured. Th interface. We as described i The addit photoelastic d undergo stretc been measure for compariso breakage [START_REF] Holmquist | The internal tensile strength of a borosilicate glass determined from laser shock experiments and computational analysis[END_REF] 4. [START_REF] Hülsenberg | Microstructuring of glasses[END_REF]

Conclusion

Femtosecond laser glass welding has been demonstrated using a long focal length lens integrated in a scanner head instead of the more conventional microscope objectives with high numerical aperture. The welding principle can be described similarly for both methods, relying on melting due to thermal accumulation effect. However, the temperature dynamics are completely different. In the case of a long focal length lens, the large focusing diameter generated a lower temperature increase by pulse. This smooth temperature dynamics reduces the residual stress induced in the welding seams and limits the risk of fracture in the glass material. The use of a long focal length lens also offers further advantages for an industrial process. The long working distance, corresponding to long Rayleigh length, reduces the difficulties involved in positioning the focusing spot on the interface. Thus, a large attainable welding scanning speed of up to 1000 mm/s has been obtained, thus being a great advantage in terms of industrialization. The mechanical and thermal resistance of the samples are besides reasonable for a welding process, with a tensile strength of 30 MPa and the ability to sustain thermal shocks of 300 °C.

  Fig

Fig. 1

 1 Fig. 1. Schemat

Fig. 3

 3 Fig. 3

Fig

  Fig. 4.

  20 beam. Due to t the laser spot ulation of this Fig. 6-a), and c ess. As the du sion inside the rival (around 1 erature of the r a microscop erature genera . 7-b).

Fig. 8 .

 8 Fig. 8. Observa

Fig. 9 .

 9 Fig. 9. Top an 5 µJ), b) top v

Fig. 10

 10 Fig. 10. Macros material abov

Table 1 .

 1

	Po
	Tole
	In most o
	focused by a
	which is alrea
	our configura
	device, it is po
	3 µJ. Increas
	welding has b
	4.2 Low resi
	The model de
	been presente

Table 1 , . Tolerance in pos ower (W) rance (µm)

 1

	of the current
	microscope ob
	ady much highe
	ation, with a h
	ossible to reach
	sing the repeti
	been obtained f
	idual stress
	eveloped to est
	ed in Section

itioning for differe 3.1 ± 75

	ereas for the F	-theta lens, the e thermal accu
	of the entire vo olume.
		0 fs, 500 kHz puls ses, showing: a) th he slow thermal ac
		h the F-theta lens,	focal length f =10 00 mm (Heat sourc
		ease with a micros scope objective, fo ocal length f =10 m
		50 μ ,	2	μ .
		phs, the use of f a long focal	length focusin
	dustrial process s: the possibili ity to reach a v
		asy positioning g of the laser b beam at the gla
		00 mm F-thet ta lens is the	long Rayleig
	Long Rayleigh h length genera ates elongated
	microscope obj jective. With s suitable scanni
	Hz, welding se eams deeper th han 200 µm h
		allows a large e freedom of	positioning ar
	with more than	±100 µm unce ertainties in po
		ocess very robu ust.
		ent laser mean po ower at a scanning
		4.4		6.0
		± 75		± 100
	t glass weldin ng demonstrat tion using fem
	bjective, the s canning speed d is limited to
	er than the firs st demonstratio on with low rep
	high repetition	rate laser syst tem and a long
	h a scanning sp peed of up to 2 200 mm/s with
	ition rate to 2 2 MHz and th he pulse energ
	for scanning sp peeds of up to 1 1000 mm/s.
	timate the evo lution of the te emperature dis
	3. Compared t to the tempera ature distributi

g speed of 50 mm 12.0 ± 200

	welding with	umulation
	induced by a l
	The residu
	in the case of
	stresses are ob
	residual stress
	where the sam
	orthogonal to
	parallel or ort
	when the seam
		cumulation
		ce:
		mm (Heat
		ng device
		very high
		ass plates
		gh length
		d welding
		ing speed
		have been
		round the
		sitioning,
		m/s
		mtosecond lase er pulses
		around 1-10 m mm/s [1],
		petition rate las ser [7]. In
		g focal length	focusing
		low energy, i.	e. around
		gy to 10 µJ, s successful
		stribution evol lution has
		ion simulated	for glass

  High qua

		microscope ob bjective, the sp peed of the tem mperature incr reases and the	gradients
		long focal leng gth focusing de evice are clearl ly much smalle er.
		ual stresses ind duced by the w welding process s [10] are thus s expected to b be smaller
	f the scanner he ead. The local b birefringence i in the glass int troduced by the e residual
		bserved by pho otoelasticimetry y, as described d in Section 2.3 3. The orientati ion of the
		ses can be det termined by th he rotation of t the sample, as s can be seen i in Fig. 8,
		me area is ob bserved, rotate d at 45°. Whe en the welded d seams are p	arallel or
	o the crossed p polarizer-analy zer system, th he field appear rs dark: the str resses are
		thogonal to the e incident pola arization. On th he opposite, th he contrast is m maximum
	ms direction is	rotated by 45	°: the stresses a are at 45° to th he incident pola arization.
		ation by photoelast ticimetry of weldin ng seams: determi ination of the resid dual stress orientat tion by the
		ro otation of the samp ple. (500 kHz, 100	mm/s, 3-6 µJ)
		tion of a firs st-order wavep plate, accordin ng to the stan ndard test me ethod for
	determination o of residual stre ess (C 978-04	standard), sho ows that weldin ng seams
		ching stress or rthogonal to th he scanning di irection [10].	The residual s stress has
	ed lower than 6 60 MPa for low w energy close e to the absorpt tion threshold	[10]. Just
		on, this is we ll below the i internal residu ual stress amou unt of 1 GPa	inducing
	].
	ality welding s seams
	The presence e of cracks a along or at th he edges of th he welding se eams perform med using
	microscope o bjective can b e induced by e either an exces ss or an overla ap of the residu ual stress
	inside the we elding seams	[1]. On the o other hand, the e low residua al stress proce ess, made
	possible by th he use of a lon ng focal length	focusing devic ce, allows for	lower risks of f fractures
	and defaults	in the welding g seams. By a adapting the p pulse energy, i it is possible	to obtain
	homogeneous s and crack-fre e welding seam ms. The side vi iew of welding g seams obtain ned by the
	scanner head	(Fig. 9) shows s a homogeneou us section with hout visible def faults.

Funding

This work has been supported by French ANRT.


Juan Carlos Quezada
email: juan-carlos.quezada@insa-strasbourg.fr

Cyrille Chazallon

Complex modulus modeling of asphalt concrete mixes using the Non-Smooth Contact Dynamics method

Keywords: Complex modulus test, asphalt concrete modeling, NSCD method, Burger's model

Asphalt mixtures are complex multiphase materials showing a viscoelastic behavior. To assess the mechanical response of these materials, a typical practice is to perform a complex modulus test. An alternative to laboratory characterization is to simulate numerically this test by means of a discrete approach. These approaches are able to reproduce the mechanical performances of asphalt mixtures, but are still time-consuming. In this paper, the complex modulus test is reproduced numerically for a viscoelastic granular material by means of 3D Non-Smooth Contact Dynamics simulations. A viscoelastic phase surrounding the rigid particles is simulated by a contact model acting between them. This contact law was implemented in the LMGC90 software, based on the Burger's model. The developed contact model handles larger time step lengths to reduce the computational time. Experimental and numerical testing campaigns were conducted for the complex modulus test on trapezoidal samples in a 2PB configuration. The numerical model was able to reproduce the mechanical performances obtained during experimental tests, regarding the material properties such as the complex modulus norm and the associate phase angle. The proposed model can be used to simulate the mechanical response of road structures under traffic loading concerning rutting, crack propagation and fatigue damage.

Introduction

Geomaterials are widely used as filling materials in construction foundations, earthworks and transport infrastructures such as railway ballast and asphalt concrete. The latter is a complex multiphase mixture comprising mastic of bitumen and filler, graded mineral aggregates and air. This mixture links the granular phase, composed of crushed aggregates, with bitumen to provide cohesion to the mixture. This material has a typical viscoelastic behavior, where the mixture properties are function of strain rate and temperature. In these materials, the macroscopic behavior is strongly linked to the micro-scale behavior. Hence, to make accurate predictions regarding pavement performances, a clear understanding of the mechanical behavior of asphalt mixtures at the particle-scale is needed.

Complex modulus test allows quantifying the viscoelastic properties of asphalt mixtures. Among these properties, the most important are the norm of the complex modulus E * together with the associate phase angle Φ. These properties are used as inputs to calculate the asphalt pavement response under traffic loading and to design the thickness of asphalt concrete pavement layers. Laboratory tests are able to identify the macro-scale response of asphalt mixtures under cyclic loading, but this approach struggles to obtain a micromechanical insight of these materials. Numerical simulations using a continuum mechanics approach such as Finite Element Method (FEM) have been successfully used to describe the overall behavior of viscoelastic materials [START_REF] Lu | Numerical approach of visco-elastoplastic analysis for asphalt mixtures[END_REF][START_REF] Collop | Development and finite element implementation of stress-dependent elastoviscoplastic constitutive model with damage for asphalt[END_REF][START_REF] Masad | Development of an elastoviscoplastic microstructural-based continuum model to predict permanent deformation in hot mix asphalt[END_REF]. Nevertheless, these methods show the same difficulties as experimental approach regarding the identification of microscale properties.

A workaround for this issue is to use a discrete approach, which has been widely used to model the behavior of granular materials. This method allows simulating the interaction of a collection of rigid or deformable bodies in contact. Over the past two decades, the discrete element method (DEM) has been used in several studies to model the mechanical behavior of asphalt mixtures. This approach takes into account the granular nature of aggregates and the mastic behavior applying viscoelastic laws at particle contact [START_REF] Mcdowell | A dimensional analysis of scaling viscosity and velocity in dem of constant strain rate tests on asphalt, Geomechanics and Geoengineering[END_REF][START_REF] Yu | A micromechanical based three-dimensional dem approach to characterize the complex modulus of asphalt mixtures[END_REF][START_REF] Dondi | Modeling the dsr complex shear modulus of asphalt binder using 3d discrete element approach[END_REF][START_REF] Feng | Study of the internal mechanical response of an asphalt mixture by 3-d discrete element modeling[END_REF][START_REF] Nguyen | Modélisation numérique discrète des matériaux bitumeux[END_REF]. These studies have modeled the asphalt using simple spherical particles, irregular particles created with clumps of spheres or by generating samples using imaging techniques [START_REF] Collop | Modelling dilation in an idealised asphalt mixture using discrete element modelling[END_REF][START_REF] Adhikari | 3d microstructural models for asphalt mixtures using x-ray computed tomography images[END_REF][START_REF] Liu | Visualization and simulation of asphalt concrete with randomly generated three-dimensional models[END_REF][START_REF] Cai | Discrete element visco-elastic modelling of a realistic graded asphalt mixture[END_REF].

Despite the accurate micromechanical behavior predictions of asphalt mixtures using classic DEM approach, this numerical method is still very time-consuming. This issue is due in part to the length of the time step dt, where this value must be smaller than a critical time step, which is determined by the mechanical properties of the contact law. In this method, increasing the time step value leads to deal with unbalanced forces, where the numerical resolution stability of the equations of motion is not ensured [START_REF] Mcdowell | A dimensional analysis of scaling viscosity and velocity in dem of constant strain rate tests on asphalt, Geomechanics and Geoengineering[END_REF].

To face these issues, in this paper a viscoelastic contact law is developed and implemented in a particle-based code as an alternative to model asphalt mixtures dynamic properties. The proposed contact model deals with larger time step values to reduce the computational time and assess the micromechanical properties of the model from macro-scale parameters at the contact particle. To validate the method and to determinate the contact model parameters, laboratory complex modulus tests in a two-point bending configuration were performed on trapezoidal samples.

The initial section (Sect. 2) of this paper set forth the experimental set-up and the determination of viscoeslatic properties from laboratory tests. Then, the numerical procedures are outlined in Sect. 3, with the description of the numerical approach, the implementation of the contact model (Sect. 3.1) and the complex modulus test simulations (Sect. 3.2). Finally, Sect. 4 analyzes the results from numerical complex modulus test and compares with experimental data to validate the proposed contact model. A sensitivity analysis of model parameters is also performed in this section.

Laboratory work

This section presents experimental procedures for the complex modulus test with the evaluation of viscoelastic macroscopic parameters.

Complex modulus set-up

Complex modulus test were conducted in a two-point bending (2PB) configuration according to the EN 12697-26:2012 specification [START_REF]Bituminous mixtures. Test methods for hot mix asphalt. Part 26: Stiffness[END_REF]. Four trapezoidal samples were prepared with Semi Coarse Asphalt Concrete (BBSG). Table 1 displays the particle size distribution (PSD) for this asphalt mixture. These aggregates come from the quarry of Bréfauchet in France and the bitumen is a 35/50 grade. The specimens were saw-cut from slabs prepared in laboratory, with the resulting dimensions: b = 25mm, B = 56mm, e = 25mm and h = 250mm (Fig. 1a). The mass of each trapezoidal sample M is on average 0.6kg and the mobile equipment mass at the head of the specimen m is 0.235kg.

For this series of trials the frequency values tested were set at 3, 6, 10, 25 and 40Hz while the selected temperatures were -10, 0, 10, 15, 20 and 30 o C. For each couple temperature-frequency, a sinusoidal displacement z = A sin(ωt) is applied on the top of each sample. The amplitude A is set to 63×10 -6 m in to obtain a strain ε less than 50×10 -6 to avoid fatigue damage. Each test is conducted for a total duration ranging from 30s to 120s. The peak values of force F 0 , displacement z 0 and initial phase angle Ψ were measured during the last ten seconds of each test.

Determination of viscoelastic material properties

The viscoelastic behavior of an asphalt mixture under cyclic loading can be assessed by the norm of the complex modulus E * and the associate phase angle Φ. To determinate these viscoelastic parameters, the real part E 1 and the imaginary part E 2 of the complex modulus must be calculated for each couple temperature-frequency. These values are provided by Eqs. 1 and 2 respectively [START_REF]Bituminous mixtures. Test methods for hot mix asphalt. Part 26: Stiffness[END_REF]:

E 1 = γ F 0 z 0 cos(Ψ) + 10 -6 µω 2 (1)
E 2 = γ F 0 z 0 sin(Ψ) (2)
γ = 12h 3 e(B -b) 3 2 - b 2B b B - 3 2 -ln b B ( 3
)
where γ (mm -1 ) is the form factor depending on the geometrical values showed in Fig. 1a and µ = 0.135M+m is the mass factor influencing the resultant force by inertial effects. Finally, the norm of the complex module and the associate phase angle are calculated using the following expressions:

E * = E 2 1 + E 2 Φ = arctan(E 2 /E 1 )
(5)

Numerical procedures

This section presents the steps to develop a contact model able to simulate a viscoelastic granular material based on the Non-Smooth Contact Dynamics (NSCD) method. Then, for the numerical modeling of the complex modulus test, the preparation protocol and the determination of model parameters are detailed.

Contact model description

In the aim to model an asphalt mixture during a complex modulus test, a viscoelastic contact model was developed and implemented in the LMGC90 software. The latter is capable of modeling a collection of rigid or deformable particles of various shapes by different algorithms [START_REF] Dubois | Actes du sixieme colloque national en calcul des structures[END_REF][START_REF] Perales | A nonsmooth contact dynamics-based multi-domain solver: code coupling (xper) and application to fracture[END_REF][START_REF] Topin | Micro-rheology of dense particulate flows: Application to immersed avalanches[END_REF][START_REF] Rafiie | The non-smooth contact dynamics method applied to the mechanical simulation of a jointed rock mass[END_REF][START_REF] Azéma | Packings of irregular polyhedral particles: Strength, structure, and effects of angularity[END_REF][START_REF] Cantor | Numerical simulation of the compaction of crushable grains in 3d[END_REF], (see https://git-xen.lmgc.univ-montp2.fr/ lmgc90/lmgc90_user/wikis/home). The created contact model mixes the original NSCD formulation [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Moreau | Some numerical methods in multibody dynamics: Application to granular materials[END_REF][START_REF] Jean | The non smooth contact dynamics method[END_REF][START_REF] Brogliato | Nonsmooth mechanics[END_REF][START_REF] Radjai | Contact dynamics as a nonsmooth discrete element method[END_REF] for rigid particles to ensure the non-interpenetration criterion with a viscoelastic model acting on distant contacts.

Non-smooth complementarity relations

For two particles candidate to the contact, it is possible to define a local frame at the particle-scale. The local frame is defined by the normal unit vector n, which links the normal gap δ n formed by the closer points between two particles. It is possible to associate to this normal vector two additional vectors of its tangential plane t and s to constitute an orthonormal base. In this local frame ( n, t, s), vector δ = (δ n , δ t , δ s ) describes the normal and tangential gaps between two particles.

In the local frame, it is possible to describe the kinematic of rigid particles at contact as complementarity relations between velocities and contact forces. For a frictional contact, these relations can be expressed by the Signorini relation (Eq. 6) and the Coulomb friction law (Eq. 7), as following:

δ n > 0 ⇒ R n = 0 δ n = 0 ⇒ U n ≥ 0, Rn ≥ 0, U n R n = 0 (6)                δ n > 0 ⇒ R t = 0 δ n = 0          U t > 0 ⇒ R t = -µR n U t = 0 ⇒ -µR n ≤ R t ≤ µR n U t < 0 ⇒ R t = µR n (7)
where U n and U t are the normal and tangential relative velocities, R n and R t are the normal and tangential local contact forces, and µ correspond to the coefficient of friction. This formulation for the Coulomb friction law is valid for both tangential directions t and s. Equation 6 expresses mutual exclusion of particles at contact, where particle-interpenetration is not allowed. Figure 2 displays a graphical representation of these relations. Besides the frictional interaction between rigid particles at contact, a viscoelastic behavior was included in the numerical modeling. For this purpose the Burger's model was chosen. This rheological model is able to reproduce the viscoelastic behavior of asphalt mixtures [START_REF] Betten | Creep mechanics[END_REF][START_REF] Liu | Determining burger's model parameters of asphalt materials using creep-recovery testing data[END_REF][START_REF] Zbiciak | Mathematical description of rheological properties of asphaltaggregate mixes[END_REF]. It describes finely the creep, relaxation and dynamic properties of asphalt mixtures and has become a commonly used contact model in DEM simulation of these materials [START_REF] Collop | Modelling dilation in an idealised asphalt mixture using discrete element modelling[END_REF][START_REF] Mcdowell | A dimensional analysis of scaling viscosity and velocity in dem of constant strain rate tests on asphalt, Geomechanics and Geoengineering[END_REF][START_REF] Yu | A micromechanical based three-dimensional dem approach to characterize the complex modulus of asphalt mixtures[END_REF][START_REF] Cai | Discrete element visco-elastic modelling of a realistic graded asphalt mixture[END_REF][START_REF] Dondi | Modeling the dsr complex shear modulus of asphalt binder using 3d discrete element approach[END_REF][START_REF] Feng | Study of the internal mechanical response of an asphalt mixture by 3-d discrete element modeling[END_REF]. The Burger's model comprises a Maxwell model combined in series with a Kelvin-Voigt model. Here, K m is the Maxwell stiffness, C m is the Maxwell viscosity, K k is the Kelvin-Voigt stiffness and C k corresponds to the Kelvin-Voigt viscosity. Figure 3 shows the schema of the Burger's model at the contact frame for the normal and tangential components. This viscoelastic model binds the particles whose gap distance δ n is below the minimum diameter d min from the numerical PSD. This contact model is applied at the barycenter of each linked particle.

n t R n R n R t R R n R n (a) (b)
Assuming a perfectly isotropic behavior, the tangential mechanical properties of the model following both tangential directions t or s are set as equals. Thus, tangential stiffness and viscosities K m t , C m t , K k t and C k t are estimated as a function of the corresponding normal components as:

K m t = K m s = K m n 2(1 + ν) (8)
C m t = C m s = C m n 2(1 + ν) (9)
K k t = K k s = K k n 2(1 + ν) (10)
C k t = C k s = C k n 2(1 + ν) ( 11
)
where ν is the Poisson's ratio. For the sake of simplicity, this ratio was set to 0.35. Regarding the formulation for the NSCD method, the relations between local velocities and contacts forces for this contact model can be expressed for the normal and tangential components by Eqs. 12 and 13:

               δ n > δ 0 |R n | ≥ R crit n ⇒ U n ≥ 0, R n = 0 |R n | < R crit n ⇒ U n = 0, R n < 0 δ n = δ 0 ⇒ U n ≥ 0, R n = 0 δ n < δ 0 ⇒ U n = 0, R n > 0 (12)                    δ t > 0 |R t | ≥ R crit t ⇒ U t ≥ 0, R t = 0 |R t | < R crit t ⇒ U t = 0, R t < 0 δ t = 0 ⇒ U t = 0, R t = 0 δ t < 0 |R t | ≥ R crit t ⇒ U t ≤ 0, R t = 0 |R t | < R crit t ⇒ U t = 0, R t > 0 (13)
where δ 0 is the initial normal gap δ n value at the initial equilibrium state, and R crit represent a yield force. The R crit value is computed as the product between the yield stress value and the projected cross-section of particles at contact. If one component of the local reaction R is larger than R crit value, the bond between particles is broken and the frictional contact conditions (Eqs. 6 and 7) are applied instead. Both normal and tangential yield force values R crit n and R crit t are set as equals. The projected cross-section at contact is obtained by πr 2 min , where r min is the minimum radius of two spheres at contact. For simplicity, the r min value was considered for the cross-section calculation; nevertheless this approach is quite conservative. Other models could be used instead, such as calculating the effective radius between two particles or determining the cross-section based on models analogous to the capillary bonds [START_REF] Soulie | Capillary cohesion and mechanical strength of polydisperse granular materials[END_REF][START_REF] Richefeu | Shear strength of unsaturated soils: Experiments, dem simulations, and micromechanical analysis[END_REF][START_REF] Scholtès | Micromechanics of granular materials with capillary effects[END_REF].

In Fig. 4 is provided a representation of the relations between local velocities and contact forces for the viscoelastic interactions. Thus, for the developed contact model, the viscoelastic Burger's model acts together with the frictional contact law. For the case of distant contacts (δ n > 0) the contact forces R are calculated using the Burger's model. When a frictional contact is produced (δ n = 0) the generated contact forces are added to the resultant force for each contact.

Non-smooth dynamics

The NSCD formulation allows to describe the equations of motion for the frictional contact law in a compact form as:

U + = U -+ W R + dt (14)
where vectors U -and U + are the contact relative velocities at times t and t + dt, respectively, and R + is the vector of local reaction at time t + dt. The Delassus local matrix W contains the values of the inverse reduced inertia, where these components are calculated as a function of the contact geometry and inertia parameters of two particles at contact [START_REF] Radjai | Contact dynamics as a nonsmooth discrete element method[END_REF]. In this form, the accelerations are replaced by velocity jumps defined by U + -U - and the equations of motion take the form of an equality between the change of momentum and the average reaction during the time step dt.

To consider the viscoelastic behavior in this contact model, it is necessary to add the contribution of the Burger's model into the general formulation of the NSCD method. In the Burger's model, the gap between particles δ is obtained by the sum of the local displacement in each section, i.e. the sum of Kelvin-Voigt section displacement δ k and the Maxwell section displacements δ mk and δ mc due to stiffness and viscosity parts respectively, given by: Considering an implicit time-integration scheme, the δ values are updated taking the values of the relative velocities U + at time t + dt generated by the viscoelastic interaction, as:

δ = δ k + δ mk + δ mc ( 15
) n t R n R n R t R R n crit (b) (a) R t crit R t crit
δ + = δ -+ dt U + (16)
where δ -and δ + correspond to the displacement at times t and t + dt. Applying the formulation of Eq. 16 for each displacement term introduced in Eq. 15 allows to provide the complementary relative velocities for the Burger's model, by:

U + k = δ + k -δ - k dt (17)
U + mk = δ + mk -δ - mk dt (18)
U + mc = δ + mc -δ - mc dt (19)
where U + k , U + mk and U + mc correspond to the complementary velocities at the end of the time step generated by the Kelvin-Voigt section, the Maxwell stiffness and the Maxwell viscos-ity, respectively. Incorporating these complementary velocity terms into the general NSCD formulation (Eq. 14), leads to the following expression :

U + = U -+ U + k + U + mk + U + mc + W R + dt (20)
Since the Burger's model is in series, the local reaction vector at contact R + (normal and tangential components) at the end of a time step dt can be expressed as:

R + = R + k = R + mk = R + mc ( 21
)
where R + k , R + mk and R + mc correspond to the local reaction vectors at time t + dt generated in the Kelvin-Voigt section, and in the Maxwell section by the stiffness and the viscosity respectively. Each local reaction can be obtained from Eqs. 22 through 24:

R + k i = -K k i δ + k i -C k i U + k i (22) R + mk i = -K m i δ + mk i (23)
R + mc i = -C m i U + mc i ( 24
)
where i stand for n, t or s, while K k i , C k i , K m i and C m i provide the normal and tangential values of stiffness and viscosities for the Kelvin-Voigt and Maxwell sections respectively.

Finally, using Eqs. 22 to 24 into Eq. 20, together with the complementarity relations generated by the viscoelastic part from Eqs 12 and 13, provides the modified NSCD formulation, given by:

U - i + K k i δ + k i C k i + R - i dtK m i + W ii + - 1 dtK m i + 1 C m i + 1 C k i R + i dt = 0 (25)
where R - i is the local reaction vector at the beginning of the time step. To calculate all the velocities and contact forces in the system, at each time step, all kinematic constraints implied by frictional contacts and viscoelasticity from Eqs. 6, 7, 12 and 13, are considered together with Eqs. 14 and 25. The timestepping schemes provided by Eqs. 14 and 25 are implicit but require explicit determination of the contact network at each time step.

This numerical problem is solved by means of an iterative procedure based on the nonlinear Gauss-Seidel method. The latter consists of solving a single contact problem with other contact forces being treated as known, and iteratively updating the reactions R + and velocities U + up to reach a convergence criterion. The iterations during a time step are stopped when the obtained R + values are stable regarding the update procedure.

For the developed contact model, the internal values are the initial normal gap δ 0 , the initial distance between barycenter of particles l 0 , the projected cross-section of particles at contact and the Kelvin-Voigt displacement δ k . At the beginning of the calculation process, the initial normal gap δ n is assessed and stocked as initial value δ 0 , together with the initial l 0 value, while the initial tangential gaps δ t and δ s are set to zero. According to Eq. 25, only the Kelvin-Voigt displacement δ k is involved in the formulation among the displacement terms presented in Eq. 15. The normal and tangential components of δ k are initially set to zero. Then, at the end of each time step, these components are updated using the Eq. 26 considering the R + current values.

δ + k i = R + k i dt + C k i δ - k i K k i dt + C k i (26)
At each time step, the yield stress criterion is checked. If the viscoelastic interaction breaks, the frictional contact law is applied instead. For our numerical simulations, a yield stress value of 2.7×10 6 Pa was fixed. This value is high enough to avoid fatigue damage during the test. This entire formulation allows implementing the developed viscoelastic contact model in the LMGC90 software.

Complex modulus test modeling

Sample preparation

The numerical samples are composed of spherical rigid particles, where each particle has a diameter ranging from 2mm to 10mm (see Tab. 2), based on the experimental PSD displayed in Tab. 1. Here, the graded curve was cut at 2mm to avoid modeling all the fines, which are considered as a part of the mastic, in the aim to reduce the number of particles in the system and consequently the computational time. Initially, for each particle-size, a specific number of spherical particles are created according to the particle size distribution. The generated spheres are disposed randomly within a cubic lattice with 0.25m length, 0.025m wide and 1.6m height. In this cubic lattice the length of each lattice element is 10mm to fit the largest particles in the sample. The initial bulk density for all particles is 2600kg.m -3 .

The preparation protocol consists in first pouring the initial spherical particles by gravity into a box composed of six plates that fit the cubic lattice described above (Fig. 5a). During this stage, the Coulomb friction law is applied for both particleparticle and particle-wall contacts, where the coefficient of friction is set to 0.2 for every contact. At the end of the deposition process, no granular segregation phenomenon was noticed (see Fig. 5a and5b). The granular layer obtained by this procedure is subjected to vibrations of small amplitude (6mm) by applying a vertical sinusoidal displacement on the top wall during 1s with a frequency of 30Hz, obtaining a final layer of 60mm height (Fig. 5b). The imposed vibration acts only on the top particle layers, in order to avoid any segregation phenomenon or Brazilian nut effect. The resulting layer is cut to fit the dimensions of the experimental specimen by removing the particles and plates whose barycenters are outside the experimental sample volume (Fig. 5c). Afterwards, the sample is turned upwards by switching the direction of gravity and placed over an immobile bottom plate, to reproduce the vertical position in the experimental setup. To simulate the mobile parts in the experimental device, a top plate is generated with 0.235kg mass and put over the sample head. The sample is fixed to both bottom and top plates using the developed contact model with a large value of yield stress to avoid any bond breakage. In this stage, the particles are reduced in a 10% of their volume in order to create gaps between them, to model the mastic surrounding the particles. This approach allows discretizing the system, concentrating the total mass of the sample within the particles. Thus, only the granular skeleton of the asphalt mixture is modeled, while the mastic phase is simulated by the action of the contact model. The bulk density of particles is increased to 3415kg.m -3 in order to match the total mass of 0.6kg of the experimental sample.

Then, the Burger's contact law is activated for all particle contacts, with the model parameters corresponding to the target temperature. The coefficient of friction between particles is set to µ = 0.7 which is a typical value used for rock crushed aggregates. Finally, by activating gravity in the simulations, the system is stabilized after 1000 time steps. Following this entire process, four numerical samples were prepared, with the same geometrical properties as experimental specimens. The total number of particles is about 12000 for all numerical specimens. Figure 5c shows a snapshot of a resulting sample.

Calibration of model parameters

This section discusses the procedure to assess the contact law parameters to reproduce experimental data from the complex modulus test. First a macromechanical model based on Burger's model parameters must be established to describe the evolution of the norm of the complex modulus E * and the phase angle Φ. In a viscoelastic material under cyclic loading, it is possible to determinate the complex compliance from the resulting strain and stress as [START_REF] Kim | Linear viscoelastic analysis of asphalt mastics[END_REF][START_REF] Tschoegl | The phenomenological theory of linear viscoelastic behavior: an introduction[END_REF][START_REF] Yu | A micromechanical based three-dimensional dem approach to characterize the complex modulus of asphalt mixtures[END_REF]:

D * = D (ω) 2 + D (ω) 2 (27)
Both the real and the imaginary part of the complex compliance, D and D respectively, can be obtained by:

D (ω) = 1 E m + E k E 2 k + ω 2 η 2 k ( 28
)
D (ω) = 1 ωη m + ωη k E 2 k + ω 2 η 2 k ( 29
)
where E m , η m , E k and η k correspond to the macroscopic Maxwell and Kelvin-Voigt stiffness and viscosities respectively. The macroscopic model parameters allow to determinate the norm of the complex modulus as E * = D * -1 and the phase angle by Φ = arctan(D /D ). Using the established analytical expressions for E * and Φ the best fit macroscopic parameters regarding experimental data can be estimated. The curve

K m n = πr 2 min l 0 E m n ( 30
)
C m n = πr 2 min l 0 η m n (31)
K k n = πr 2 min l 0 E k n (32)
C k n = πr 2 min l 0 η k n (33)
where r min correspond to the minimum radius of two particles in contact and l 0 is the initial distance between the barycenters of two particles (see Fig. 3). Then, tangential components can be assessed from normal ones using Eq. 8 through Eq. 11. This approach creates a number of parameters at the micro-scale equal to the number of particles in contact within the sample. Table 3 displays the set of calibrated parameters determined with this procedure for each temperature.

To ensure numerical stability in the resolution, the time step length must be smaller than the characteristic time scale of the system. First, it is possible to quantify the typical time of rearrangement of grains dt g = m g /(pd g ), with the average values of particle mass m g 2.86×10 -5 kg, pressure p 2.79×10 3 Pa, and particle size d g 2.52×10 -3 m. In our system, the dt g value is about 2×10 -3 s. On the other hand, due to the elasticity in the proposed contact model, a second time scale can be identify as dt e = m g /K eq , where K eq is the equivalent contact stiffness. In the proposed model, at the beginning of the calculation time, K eq K m . For the case of 15 o C, a dt e value of 2×10 -5 s is obtained.

The influence of the time step dt on the total reaction measured on the top plate was assessed. Figure 6 displays the curves for the total reaction as a function of time for several values of dt varying from 2×10 -5 s to 2×10 -4 s. These results show that force values can reach a stable measure with a time step value of 5×10 -5 s, even if the latter is larger than the identified time scale. The time step was set to 5×10 -5 s in all simulations and the number of time steps needed to perform five loading cycles ranging from 3.33×10 4 to 2.5×10 3 for 3Hz and 40Hz respectively. Using this time step value implies that the total computational time in these numerical simulations is divided by 2.5, compared to the typical time scale used in classic DEM approaches. The CPU time was 2×10 -4 s per particle and per time step on a Dell computer of speed 2.1GHz.

Validation of the numerical modeling

Comparison between experimental and numerical data

The numerical complex modulus tests were conducted by subjecting the top plate to the same sinusoidal displacement as experimental tests. Then, five cycles were applied for each sample, measuring the displacement and the total reaction on the top plate as a function of time. The obtained force and displacement curves allow to determinate the peaks values F 0 and z 0 and calculate E * and Φ from numerical trials. Figure 7 displays an example of the imposed displacement and the corresponding force measured at the top plate. As expected, the resultant curve of the total reaction follows the sinusoidal wave imposed by the displacement. The phase angle between these two curves allows to calculate Φ. For all samples, six temperatures and five frequencies were tested such as for experimental trials. Figure 8 depicts different stages during the numerical test. The mean velocity of particles increases as a function of the height, from zero for the particles bonded to the bottom plate, up to a maximum value for the particles in contact with the top plate. The numerical sample follows the movement imposed by the top plate as an viscoelastic beam without a local effect generated by the loading plate.

To validate the numerical model, Fig. 9 and 10 display the experimental and numerical results for the norm of complex modulus and phase angle isotherms respectively. One can note that the average numerical values from the simulations are in a good agreement with the values for E * found through experiments for all the tested frequencies. On the other hand, the numerical results for the phase angle Φ are just in a fairly good agreement with the experimental results. These isotherms show some fluctuations around the average experimental data.

In the aim to assess the forecasting error, Tab. 4 gives the results of the mean relative error (MRE) between experimental and numerical data for E * and Φ for each temperature. As we can see in Fig. 9, the numerical values of E * are very close to experimental ones, with a MRE smaller than 4% for all temperatures. On the other hand, the MRE values for Φ remain below 10% except for the case of -10 o C. These variations were already pointed out in Fig. 10. This high variability at low temperature is produced in part for a ratio of two small numbers in the calculation of MRE, both the numerical and experimental Φ values. Despite the variations in the prediction for some temperatures

Sensitivity analysis of model parameters

To complete the study and to take into account a lack of accuracy into the calibration of Burger's model parameters, Fig. 11 to 14 display the influence of these parameters in the determination of E * and Φ. This sensitivity analysis was performed for 20 o C, where each model parameter was fixed at ±50% regarding the average values given in Tab. 3.

Regarding the effect of the model parameters variation on E * values, E m shows the most important influence compared to the other parameters. Theoretically, the stiffness of the Maxwell model can be express as E m = lim f →∞ E * . Hence, the influence of E m values on E * increase as the frequency increases, as displayed in Fig. 11. A similar trend but with a lower magnitude can be highlighted for the Kelvin-Voigt viscosity η k . Concerning the influence on the phase angle, a variation of 50% on each parameter of the contact model, leads to significant fluctuations in the obtained Φ values. It seems difficult to figure out the origin of this high sensibility regarding the phase angle. One can point out the need of a correct choice of the model parameters to avoid biases in the prediction of this material property.

Concluding remarks

In this paper NSCD simulations of the complex modulus test in a 2PB configuration were performed. To model the viscoelastic phase surrounding rigid spherical particles in an asphalt mixture, a viscoelastic contact law based on the Burger's model was implemented in the LMGC90 software. These simulations reproduce the experimental campaign, where four trapezoidal samples were subjected to five frequencies and six temperatures.

The outcomes in this study show that the implemented viscoelastic contact law based on the Burger's model is able to reproduce the complex modulus properties of asphalt mixtures. For different analysis conducted in this paper, the numerical results of these simulations were in a good agreement regarding experimental data. For isotherm curves, experimental and numerical data of the norm of the complex modulus E * collapse in the same curve. The numerical phase angle Φ values follow the trend of experimental ones, despite some fluctuations around the average values. Sensitivity analysis shows that the Maxwell stiffness has an important influence in the calculation of E * , increasing with the frequency values. For the case of Φ, the contact model is quite sensitive to changes in parameter values, diverging from the average tendency.

In conclusion, this work can be considered as encouraging for the validation of the proposed contact model for the study of asphalt mixtures features such as creep, rutting and fatigue damage. In this study, the parameters of the Burger's model were calibrated using analytical expressions to assess the micromechanical properties for the contact model from macroscale parameters obtained by experimental tests.

In this preliminary work, the numerical test was modeled using spherical aggregates. In a granular material, as an asphalt mixture, the particle shape plays a major role regarding fabric anisotropy, force transmission and friction mobilization. This aspect must be taken into account when extrapolating these results to actual systems.

Further works will focus on the characterization of the rheological properties of mastic to incorporate them directly into the contact model. This feature allows forecasting the viscoelastic properties of an asphalt mixture from its formulation. The next steps in this work will be to study the mechanical behavior of multi-layer road structures and the influence of particle shape and different particle size distributions on the mechanical performances taking into account actual aggregates from 3D scan to obtain more realistic asphalt mixtures. The numerical results showed in this study were obtained for small amplitude vibrations. In the case of large deformation or failure analysis, the yield stress value must be calibrated to extend the developed model to study cases such as crack propagation or fatigue damage.

Figure 1 :

 1 Figure 1: (a) Dimensions of the trapezoidal sample. (b) Schema of the complex test in a two-point bending configuration.

Figure 2 :

 2 Figure 2: Relations between velocities and contact forces for the frictional contact law: (a) The Signorini relation (normal); (b) Coulomb friction law (tangential).

Figure 3 :

 3 Figure 3: Schema representation of the Burgers model in the local frame: (a) Normal components; (b) Tangential components.

Figure 4 :

 4 Figure 4: Relations between velocities and contact forces for the Burger model: (a) Normal components; (b) Tangential components.

Figure 5 :

 5 Figure 5: Snapshots for several stages during the preparation process. (a) Initial deposit by gravity. (b) Numerical sample before compaction. (c) Numerical sample at the end of the preparation.

Figure 6 :

 6 Figure 6: Influence of time step on the measured force on the top plate during one cycle. These simulations were performed for 15 o C and 3Hz.

Figure 7 :

 7 Figure 7: Applied displacement and corresponding force at last cycle for 15 o C and 40Hz.

Figure 8 :Figure 9 :

 89 Figure 8: Average velocity of particles as a function of time during a half of a cycle. This simulation was conducted for 20 o C and 3Hz.

Figure 10 :

 10 Figure 10: Phase angle Φ isotherms for experimental and numerical tests.

  .02e+7 [Pa] E m =1.01e+7 [Pa] E m =3.03e+7 [Pa]

Figure 11 :

 11 Figure 11: Influence of E m on the obtained values of the norm of complex modulus E * (a) and phase angle Φ (b)

Figure 12 :

 12 Figure 12: Influence of η m on the obtained values of the norm of complex modulus E * (a) and phase angle Φ (b)

  .72e+7 [Pa] E k =1.36e+7 [Pa] E k =4.08e+7 [Pa]

Figure 13 :

 13 Figure 13: Influence of E k on the obtained values of the norm of complex modulus E * (a) and phase angle Φ (b)

  .61e+5 [Pa.s] η k =2.31e+5 [Pa.s] η k =6.92e+5 [Pa.s]

Figure 14 :

 14 Figure 14: Influence of η k on the obtained values of the norm of complex modulus E * (a) and phase angle Φ (b)

Table 1 :

 1 Particle size distribution for the BBSG 0/10 Bréfauchet asphalt

	diameter (mm) PSD (%)
	12.5	100
	10	90
	8	75
	6.3	58
	4	44
	2	31
	1	20
	0.5	15
	0.25	12
	0.125	9
	0.063	6.5

Table 2 :

 2 Size and number of generated particles within the numerical samples.

	Particle size (mm)	2	4 6.3	8 10
	Number of particles 11798 1588 492 211 71

Table 3 :

 3 Burger's model parameters used in numerical simulations. fitting was performed by means of an optimization algorithm from Python programming language, for each temperature. Finally the microscopic normal model parameters (K m n , C m n , K k n and C k n ) are obtained from the macroscopic model parameters (E m n , η m n , E k n and η k n ) through:

	-10	3.98×10 7 6.28×10 7 9.56×10 8 1.49×10 7
	0	3.55×10 7 2.26×10 7 2.95×10 8 3.46×10 6
	10	2.78×10 7 6.26×10 6 8.46×10 7 9.94×10 5
	15	2.43×10 7 3.71×10 6 4.84×10 7 7.73×10 5
	20	2.02×10 7 1.79×10 6 2.72×10 7 4.61×10 5
	30	1.30×10 7 3.45×10 5 9.31×10 6 1.31×10 5

T ( o C) E m (Pa) η m (Pa.s) E k (Pa) η k (Pa.s)

Table 4 :

 4 Measure of mean relative error (MRE) in forecasting values of E * and Φ.

	T ( o C)	-10	0	10	15	20	30
	MRE in E * (%) 0.56 1.33 3.82 3.31 2.86 3.72
	MRE in Φ (%)	27.7 9.75 9.81 5.47 4.55 7.01
	regarding Φ values, the general results can be considered as en-
	couraging for a validation of the proposed contact model for the
	study of asphalt mixtures.

[START_REF] Mcdowell | A dimensional analysis of scaling viscosity and velocity in dem of constant strain rate tests on asphalt, Geomechanics and Geoengineering[END_REF]

Acknowledgments

The work presented in this paper was financially supported by the French National Research Agency (ANR -SolDuGri project ANR-14-CE22-0019). The authors want also to acknowledge CST Colas France for providing the experimental data.


Guoliang Luo
email: luoguoliang@ecjtu.edu.cn

Zhigang Deng
email: zdeng4@uh.edu

Xiaogang Jin
email: jin@cad.zju.edu.cn

Xin Zhao

Wei Zeng

Wenqiang Xie

Hyewon Seo

3D Mesh Animation Compression based on Adaptive Spatio-temporal Segmentation

Keywords: Computing methodologies → Computer graphics, Animation, 3D mesh animation, compression, adaptive spatio-temporal segmentation

With the recent advances of data acquisition techniques, the compression of various 3D mesh animation data has become an important topic in computer graphics community. In this paper, we present a new spatio-temporal segmentation-based approach for the compression of 3D mesh animations. Given an input mesh sequence, we first compute an initial temporal cut to obtain a small subsequence by detecting the temporal boundary of dynamic behavior. Then, we apply a two-stage vertex clustering on the resulting subsequence to classify the vertices into groups with optimal intraaffinities. After that, we design a temporal segmentation step based on the variations of the principle components within each vertex group prior to performing a PCA-based compression. Our approach can adaptively determine the temporal and spatial segmentation boundaries in order to exploit both temporal and spatial redundancies. We have conducted many experiments on different types of 3D mesh animations with various segmentation configurations. Our comparative studies show the competitive performance of our approach for the compression of 3D mesh animations.

INTRODUCTION

The key information of a 3D mesh animation is its dynamic behavior, which drives the deformations of different mesh surface areas. As reported in existing literature, we can achieve a better performance on the compression of 3D mesh animations with repetitive motions or rigid mesh segments, which contain significant redundancies either temporally or spatially [START_REF] Aris | Adaptive compression of animated meshes by exploiting orthogonal iterations[END_REF][START_REF] Nikolče | SPC: fast and efficient scalable predictive coding of animated meshes[END_REF][START_REF] Váša | Compressing dynamic meshes with geometric laplacians[END_REF]]. Therefore, it is important to exploit the dynamic behaviors based on both spatial and temporal segmentations within a 3D mesh animation for effective data compression. However, due to the high complexity and the large data size, it remains a challenge to jointly explore the spatial and temporal segmentations to further improve the performance of 3D mesh animation compression.

In this paper, we propose an adaptive spatio-temporal segmentation based model for the compression of 3D mesh animations. Specifically, we first introduce a temporal segmentation scheme that explores the temporal redundancy by automatically determining the optimal temporal boundaries. Then, we also introduce a novel two-stages vertex clustering approach to explore the spatial redundancy by automatically determining the number of the vertex groups with optimal intra-affinities. As an application of the above adaptive spatio-temporal segmentation model, we develop a full scheme for the compression of 3D mesh animations (Figure 1). Through many experiments, we show the effectiveness and efficiency of our approach, compared to the state of the art 3D mesh animation compression algorithms.

The contributions of this work can be summarized as follows:

• an adaptive spatio-temporal segmentation approach which explores both the spatial and temporal redundancies for 3D mesh animations; and

• a compression model for 3D mesh animations by coupling the novel adaptive spatio-temporal segmentation and the compression of 3D mesh animations. The remainder of this paper is organized as follows. We first review previous and related works on the compression of 3D mesh animations in Section 2. In Section 3, we briefly give the overview of our compression scheme. Then, we present the details of our spatio-temporal segmentation model and its application to the compression of 3D mesh animations in Section 4. The experimental results by our model are shown in Section 5. Finally, we conclude this work in Section 6.

RELATED WORK

The compression of 3D mesh animation data has been a persistent research topic in the past several decades [START_REF] Maglo | 3D Mesh Compression: Survey, Comparisons, and Emerging Trends[END_REF]].Among the existing methods, a large portion of the methods take a matrix form of the 3D mesh animation, on which many of classical data compression methods and algorithms can be applied, including Principal Component Analysis (PCA) [START_REF] Alexa | Representing Animations by Principal Components[END_REF][START_REF] Hou | Sparse Low-Rank Matrix Approximation for Data Compression[END_REF][START_REF] Liu | Limited Memory Block Krylov Subspace Optimization for Computing Dominant Singular Value Decompositions[END_REF], linear prediction encoders [START_REF] Karni | Compression of soft-body animation sequences[END_REF][START_REF] Stefanoski | Scalable Linear Predictive Coding of Time-Consistent 3D Mesh Sequences[END_REF][START_REF] Nikolče | SPC: fast and efficient scalable predictive coding of animated meshes[END_REF][START_REF] Hyu | Compression of 3-D triangle mesh sequences based on vertex-wise motion vector prediction[END_REF], wavelet decomposition [START_REF] Guskov | Wavelet Compression of Parametrically Coherent Mesh Sequences[END_REF][START_REF] Payan | Temporal wavelet-based compression for 3D animated models[END_REF], and the Moving Picture Experts Group (MPEG) framework [START_REF] Mamou | Framebased compression of animated meshes in MPEG-4[END_REF]. PCA is a classical method that can decompose a large matrix as the product of two much smaller matrices, with minimal information loss. Following the work of [START_REF] Alexa | Representing Animations by Principal Components[END_REF], [START_REF] Lee | 3D animation compression using affine transformation matrix and principal component analysis[END_REF] apply PCA to 3D mesh animation data after removing its rigid transformations. Later, researchers have used the linear prediction theory to further encode the resulting coefficients from PCA [START_REF] Karni | Compression of soft-body animation sequences[END_REF][START_REF] Vasa | COBRA: Compression of the Basis for PCA Represented Animations[END_REF][START_REF] Váša | Exploiting Connectivity to Improve the Tangential Part of Geometry Prediction[END_REF]. Similarly, researchers have proposed a Laplacian-based spatio-temporal predictor [START_REF] Váša | Compressing dynamic meshes with geometric laplacians[END_REF] or curvature-and-torsion based analysis [START_REF] Yang | Motion-aware Compression and Transmission of Mesh Animation Sequences[END_REF] to encode the vertex trajectories for dynamic meshes. However, they assume an entire sequence as the given input, and do not explicitly exploit the dynamic behaviors enclosed in the input animation. The key information of a 3D mesh animation is its enclosed dynamic behavior; therefore, it is important to exploit the dynamic behavior coherence in 3D mesh animations for effective compression, using either spatial segmentation or temporal segmentation methods.

Spatial segmentation based compression:

The key of the spatial segmentation of a 3D mesh animation is to understand its semantic behaviors. Many previous methods have been proposed to compute the spatial segmentations for 3D mesh animations, which can generate different spatial segmentation schemes for animations with different motions [START_REF] Vasilakis | Pose partitioning for multi-resolution segmentation of arbitrary mesh animations[END_REF][START_REF] Edilson De Aguiar | Automatic Conversion of Mesh Animations into Skeleton-based Animations[END_REF][START_REF] James | Skinning Mesh Animations[END_REF][START_REF] Kavan | Fast and Efficient Skinning of Animated Meshes[END_REF][START_REF] Huy | Robust and accurate skeletal rigging from mesh sequences[END_REF][START_REF] Tong | Segmenting a deforming mesh into near-rigid components[END_REF][START_REF] Wuhrer | Segmenting animated objects into near-rigid components[END_REF]. [START_REF] Hijiri | A spatial hierarchical compression method for 3D streaming animation[END_REF] separately compress the vertices of each object with the same movements to obtain an overall optimal compression rate. In order to adapt spatial segmentation for compression, [START_REF] Sattler | Simple and efficient compression of animation sequences[END_REF] proposed an iterative clustered PCA method to group the vertex trajectories that share similar Principal Component (PC) coefficients and then further compress each cluster separately. Its main limitation is its heavy computational cost. Similarly, [START_REF] Subramanian Ramanathan | Impact of vertex clustering on registration-based 3D dynamic mesh coding[END_REF] compute the optimal vertex clustering for the optimal compression ratio. However, all the above methods assume the entire animation has been given at the beginning.

Temporal segmentation based compression:

The objective of temporal segmentation is mainly to chop a 3D mesh animation into sub-sequences, each of which represents a different dynamic behavior. Temporal segmentation has been exploited for the compression of motion capture data [START_REF] Gong | Kernelized temporal cut for online temporal segmentation and recognition[END_REF][START_REF] Gu | Compression of human motion capture data using motion pattern indexing[END_REF][START_REF] Sattler | Simple and efficient compression of animation sequences[END_REF][START_REF] Zhu | Quaternion space sparse decomposition for motion compression and retrieval[END_REF]], but the efficiencies of these methods for 3D mesh animation compression may be significantly decreased since 3D mesh surfaces typically have much more denser vertices and additional topology than motion capture data [START_REF] Luo | Joint entropy-based motion segmentation for 3D animations[END_REF]. Given a mesh sequence, [START_REF] Luo | Compression of 3D mesh sequences by temporal segmentation[END_REF] group the meshes with similar poses and apply PCA to compress each group to achieve the optimal compression ratio. Recently, Lalo et al. [2017] proposed an adaptive Singular Value Decomposition (SVD) coefficient method for 3D mesh animation compression. They first divide a mesh sequence into temporal blocks of the same length and treat the first block with SVD. Then, the following blocks are treated with the adaptive bases from the previous block without solving the full SVD decomposition for each block, which reduces the compute time.

In summary, spatial and temporal segmentations can help to reveal the spatial and temporal redundancies within 3D mesh animations, which benefits for the development of effective compression algorithms. The new compression scheme for 3D mesh animations, presented in this work alternately exploits both spatial and temporal redundancies.

SCHEME OVERVIEW

In general 3D mesh animations mainly have two different forms, namely, time-varying meshes and deforming meshes. A time-varying mesh may have different numbers of vertices and different topological connectivities at different frames, whereas a deforming mesh has a fixed topology across frames. Note that we can always compute the inter-frame vertex correspondences to convert a time-varying mesh into a deforming mesh [START_REF] Tevs | Animation Cartography -Intrinsic Reconstruction of Shape and Motion[END_REF]. For the sake of simplicity, we focus on the deforming mesh data in this work.

Then, we define the trigger conditions for the two important steps in our method. (1) Initial Temporal Cut: given the maximal length γ init if any dynamic behavior has been detected in the mesh sequence (with no more than γ init frames) (see Section 4.1), and (2) Actual Temporal Segmentation: given the maximal length γ act if any dynamic behavior has been detected in any of the vertex groups (see Section 4.2 and 4.3).

We briefly describe the pipeline of our segmentation scheme as follows. The algorithmic description is also shown in Figure 1.

(1) We first conduct an initial temporal cut to produce a subsequence S with the maximal possible length of γ init , see Section 4.1.

(2) If no distinct behavior can be detected in S, i.e., the boundary frame b = γ init , the subsequence S will be directly sent to the compressor (the case (I) in Section 4.5), see Section 4.4. (3) Otherwise (i.e., distinct behaviors are detected in S), we perform a 2-stages vertex clustering on S, see Section 4.2. (4) Then, we continue to compute the temporal segmentation of each vertex group (spatial segment) within next γ act frames, by analyzing the dynamic behaviors, see Section 4.3. (5) If we have detected distinct dynamic behaviors of any vertex group before γ act is reached, the vertex trajectories of each group up to the detected boundary frame are sent to the compressor, separately. See Section 4.4. After the compression, we repeat the process from the step 1 (the case (II) in Section 4.5). ( 6) Otherwise (i.e., we have not detected a temporal segmentation before reaching γ act ), we also send the data of each vertex cluster to the compressor, separately (the case (III) in Section 4.5). See Section 4.4. Afterwards, we reuse the previously obtained vertex clustering and continue the analysis of the temporal segmentation in the remaining mesh frames. That is, we repeat the process from the step 4 for the remaining mesh frames (The case (IV) in Section 4.5).

SPATIO-TEMPORAL SEGMENTATION FOR COMPRESSION

We first describe our spatio-temporal segmentation model that consists of the initial temporal cut (Section 4.1), vertex clustering (Section 4.2), and temporal segmentation (Section 4.3). Then, we apply the spatio-temporal segmentation model for the compression of 3D mesh animations in Section 4.4. Finally, we discuss different scenarios while processing a continuous mesh sequence as the input in Section 4.5.

Initial Temporal Cut

Let us denote a mesh animation as ({V

f i }, E)
, where E represents the connectivities among the vertices, and

V f i = (x f i , y f i , z f i ) represents the 3D coordinates of the i-th vertex (i = 1, . . . , V ) at the f -th frame (f = 1, . . . , F ).
Here V is the total number of vertices, and F is the total number of frames in the animation sequence.

Given a mesh sequence, the objective of the initial temporal cut is to determine a boundary frame

V |τ | , so that the dynamic behavior in [V 1 , V |τ | ] is distinctive from that in [V |τ |+1 , V γ ini t ].
To this end, we can formulate the initial temporal cut as the following optimization problem:

min b ∈[1,γ ini t ] I ([V 1 , V b ], [V b+1 , V γ ini t ]), (1)
where b is a to-be-solved frame index and I (•, •) computes the affinity between two mesh subsequences. Available techniques for computing I (•, •) can be classified into two categories: 1) front-to-end, uni-directional boundary candidate search, and 2) bi-directional boundary candidate search. Between them, the bi-direction search method is more robust on detecting the temporal cut between two successive dynamic behaviors [START_REF] Barbič | Segmenting Motion Capture Data into Distinct Behaviors[END_REF][START_REF] Gong | Kernelized temporal cut for online temporal segmentation and recognition[END_REF]. Inspired by the kernelized Canonical Correlation Analysis (kCCA) approach [START_REF] Hofmann | Kernel methods in machine learning[END_REF][START_REF] Smola | A Hilbert space embedding for distributions[END_REF], and its successful application to semantic temporal cut for motion capture data [START_REF] Gong | Kernelized temporal cut for online temporal segmentation and recognition[END_REF], we formulate the initial temporal cut to a Maximum-Mean Discrepancy problem as follows:

min b i ∈[1,γ ini t -ϵ ] - 1 |T 1 | 2 |T 1 | i, j K(v b i :b i +ϵ , v b j :b j +ϵ ) -1 |T 1 | |T 2 | |T 1 | i |T 2 | j K(v b i :b i +ϵ , v b j :b j +ϵ ) + 1 |T 2 | 2 |T 2 | i, j K(v b i :b i +ϵ , v b j :b j +ϵ ) , (2) whereT 1 is the subsequence [V 1 , . . . , V b i ] andT 2 is the subsequence [V b i +1 , . . . , V γ ini t -ϵ ],
and ϵ is a pre-defined parameter to ensure smooth kernels.

The kernel function in Eq. 2 is defined as follows:

K(v b i :b i +ϵ , v b j :b j +ϵ ) = exp(-λ∥v b i :b i +ϵ -v b j :b j +ϵ ∥ 2 ), ( 3
)
where λ is the kernel parameter for K(•) [START_REF] Van Vaerenbergh | Kernel methods for nonlinear identification, equalization and separation of signals[END_REF].

Due to the symmetric property of the kCCA, i.e., K(A, B) = K(B, A), we obtain a symmetric kCCA matrix for the animation block. Finally, we can obtain a boundary frame V |τ | for the initial temporal cut by solving the objective function in (Eq. 2). Note that |τ | = b + ϵ due to the usage of a smoothing window. Meanwhile, we denote the detected initial temporal cut as τ . Figure 2 shows one of the initial temporal cuts of the 'March' data, with γ init = 20 and ϵ = 5. The complexity of the above bi-directional search for the initial temporal cut is O(|γ init | 2 ), which is less efficient than the unidirectional methods with O(|γ init |). However, in our context, we compute the initial temporal cut within a short mesh sequence [1, γ init ], which is a small cost on the computation and thus will not cause notable delay to the overall compression framework. The settings of γ init for different experimental data are presented in Table 1.

Vertex Clustering

In this section, we describe a vertex clustering (spatial segmentation) algorithm based on a two-stages, bottom-up hierarchical clustering algorithm to obtain optimal spatial affinities within segments.

Initial Vertex

Clustering. After the initial temporal cut, τ is obtained; we then compute the vertex clustering based on the dynamic behaviors of different vertices. The pipeline of our approach is shown in Figure 3 (I,II,III).

In this initial vertex clustering stage, we first segment a dynamic mesh based on rigidity with the following steps.

(1) Compute the MEC for all edge pairs. Similar to [START_REF] Tong | Segmenting a deforming mesh into near-rigid components[END_REF] Then, with the aid of the inverse cumulative distribution function of epd, we can determine a user-specified percent of the edges as the rigid edges (ρ = 20% in our experiments). Thus, the vertices that are connected to the rigid edges are Inspired by this work, we design the second-stage vertex clustering by iteratively classifying and assigning each rigid cluster to the group with the minimal reconstruction error until the grouping remains unchanged. Since the iterative clustered-PCA is performed on the initial vertex cluster, it works very efficiently, unlike the case in [START_REF] Sattler | Simple and efficient compression of animation sequences[END_REF]].

( Ⅰ ) ( Ⅱ ) ( Ⅲ ) ( Ⅳ ) 'deformed' 'rigid'
The reconstruction error of a rigid cluster δ j is defined as follows:

∥δ j -δ j ∥ 2 = ∥δ j -(C[j] + δ j )∥ 2 , ( 4
)
where δ j is the reconstructed cluster using PCA, C[j] is the center of each group (j = 1, . . . , ω), and δ j is the reconstruction using the PCA components (see Eq. 6). Note that we have C[j] in Eq. 4, because PCA contains the centering (mean subtraction) of the input data for the covariance matrix calculation. Figure 3 illustrates the process of the rigid cluster grouping with the 'March' data. As an example result in Figure 3(IV), the relatively less moved rigid clusters, 'head', 'chest' and 'right-arm', are classified into the same group. Note that we obtain large vertex groups because the input mesh are smooth on the surface (see Table 1), unlike the motion Capture data containing sparse vertex trajectories that may lead to small groups. Moreover, the computational cost for the initial vertex clustering presented above is relatively small because both the number of clusters k and the number of groups ω are small.

Temporal Segmentation

After obtaining a set of the spatio-temporal segments L(δ ) j (j = 1, . . . , ω) for the initial temporal cut τ , we further introduce a temporal segmentation step as follows:

• For each vertex group, we stop observing the number of PCs once it is changed within the current sliding window. In this way, we can obtain a Num-of-PCs curve for each vertex group, see the bottom of Figure 4. • To this end, similar to [START_REF] Karni | Compression of soft-body animation sequences[END_REF], the temporal segmentation boundary is determined as the first frame where any Num-of-PCs curve has changes, see the bottomright of Figure 4.

The complexity of the temporal segmentation is O(ωγ act ), where γ act denotes the maximal length of temporal segments. Note that the computational cost of the PCA decomposition increases exponentially with the input data size. In order to balance the computational cost and the effectiveness of PCA, we set an adaptive γ act for each of the input data, see Table 1. Parallel computing. The temporal segmentation presented above is designed for each vertex group (spatio-temporal segment), and the vertex groups are independent of each other. Thus, we can implement the temporal segmentation for each vertex group in parallel. The computational time statistics in Table 1 show the efficiency improvement through parallelization.

Compression

After the above spatio-temporal segmentation, we apply PCA to compress each segment with a pre-defined threshold on the information persistence rate, µ ∈ [0 1], which is used to determine the number of PCs to retain after the PCA decomposition, i.e.,

k i (σ i )/ |n | i (σ i ) ≥ µ, (5)
where k ≤ n, and {σ i }(i = 1, . . . , n) are the eigen-values of the data block in a descending order. Therefore, we can control the compression quality by manipulating the value of µ. In specific, by increasing µ, we have less information loss but more storage costs after compression; and vice-versa.

• Encoder. For a spatio-temporal segment L(δ ) i j , i.e., the j-th spatial segment within the i-th actual temporal segment ατ i , we denote its compression as follows:

V α τ i L(δ ) i j PCA ≈ A i j × B i j , (6)
where A i j is the score matrix of dimensions 3|V L(δ

) i j | × k i j , B i j is the coefficient matrix of dimensions k i j × |ατ i |,
and X denotes a centered matrix of X by subtracting the mean vectors X, i.e., X = X -X.

(7)

• Decoder. With the score matrix and the coefficient matrix, we can approximate each of the spatio-temporal segments using Eq. 6 and Eq. 7. Then, we can reconstruct the original animation by concatenating the spatio-temporal segments in order.

Sequential Processing

As discussed in Section 3, our spatio-temporal segmentation scheme generates four possible animation blocks that are further sent to the encoder for compression (see Figure 1), which leads to four types of sequential processing to the successive mesh sequence:

(I) |τ | = γ init . This indicates none of distinct behaviors has been detected at the initial temporal cut step (Section 4.1). In this case, the animation block [V 1 , V γ ini t ] will be directly sent to the encoder. Moreover, we need to re-compute a spatio-temporal segmentation for the successive mesh sequence.

(II) |τ | < γ init and |ατ | < γ act . This indicates the vertex clustering has been conducted and a temporal segmentation boundary has been detected at V |α τ | . In this case, each vertex group of the animation block will be sent to the encoder, separately. Moreover, we will re-compute a spatio-temporal segmentation for the successive mesh sequence.

(III) |τ | < γ init and |ατ | = γ act . This indicates the vertex clustering has been conducted and a temporal segmentation boundary has not been detected within the range [V 1 , V γ ac t ] . In this case, each vertex group of the animation block will be sent to the encoder, separately. Moreover, we will only need to re-compute the temporal segmentation for the successive mesh sequence.

(IV) Otherwise, we can directly reuse the existing (previous) vertex grouping results, compute the temporal segmentation, and then perform the PCA-based compression for each vertex group. If the new boundary |ατ | < γ act , we will need to re-compute a Table 1: The results and performances by our model with different configurations of parameters: ϵ and γ init for the Initial Temporal Segmentation (Section 4.1), γ act for the Actual Temporal Segmentation (Section 4.3) and ω for the vertex clustering (Section 4.2). s and s p are the timings in seconds (unit) of the single-thread and paralleled implementations, respectively, with the last column showing the percentage of the time savings for each data. spatio-temporal segmentation for the successive mesh sequence; otherwise (i.e., |ατ | = γ act ), it will again become the case (IV) for the successive mesh sequence.

Animations

EXPERIMENT RESULTS AND ANALYSIS

In this section, we first present the experimental data and the used evaluation metrics in Section 5.1. Then, we describe our experimental results in Section 5.2. In addition, we conducted a comparative study in Section 5.3. Both our model and the comparative methods were implemented with Matlab and the experiments were performed on an Intel Core i5-6500 CPU @3.2GHz (4 cells) with 12G RAM. More results can be found in the supplemental materials.

Experimental Setup

Table 1 shows the details of our experimental data. Among them, 'March', 'Jump' and 'Handstand' were created by driving a 3D template with multi-view video [START_REF] Vlasic | Articulated Mesh Animation from Multi-view Silhouettes[END_REF]]. 'Horse' was generated by deformation transfer [START_REF] Sumner | Deformation Transfer for Triangle Meshes[END_REF]. 'Flaд' and 'Cloth' are dynamic open-edge meshes [START_REF] Cordier | A Data-Driven Approach for Real-Time Clothes Simulation[END_REF]. We applied the following two metrics for quantitative analysis: Bits per vertex per frame (bpvf). Similar to [START_REF] Chen | Cloth compression using local cylindrical coordinates[END_REF][START_REF] Nikolče | SPC: fast and efficient scalable predictive coding of animated meshes[END_REF], we also used bpvf to measure the performance of compression methods. By assuming that the vertex coordinates are recorded as single-precision floating numbers, the bpvf of the original animation is 8bits/Byte × 4Bytes × 3 = 96. Thus, we can estimate the bpvf of our model as follows:

bpv f = 96 • i, j (3 • |V L(δ ) i j | × k i j + k i j × |ατ i | + |ατ i |)/(3 • V × F ).
(8) Reconstruction error. After compression, we can reconstruct the animation with the decoder described in Section 4.4. In order to measure the difference between the reconstructed animation and the original animation, we use the well-known metric KGError, proposed by [START_REF] Karni | Compression of soft-body animation sequences[END_REF]:

100 • ∥F -F∥ f ∥F -E(F)∥ f , (9)
where ∥ • ∥ f denotes the Frobenius norm, F and F are the original animation coordinates and the reconstructed animation coordinates of size 3V × F , respectively. Moreover, E(F) are the averaged centers of each frame, and thus F -E(F) indicates the centering of the original animation.

Experimental Results

We present and discuss both the segmentation results and the compression results in this section.

Spatio-temporal segmentation results. Figure 5 shows some samples of the spatio-temporal segmentation results of our experimental data (more results can be found in our supplemental materials). As can be seen in this figure, given the maximal number of spatial segments (groups) ω = 4, our model is able to automatically determine the optimal number of vertex groups (i.e., exploiting the spatial redundancy) for different dynamic behaviors (i.e., exploiting the temporal redundancy) for all the data. For example:

• The segmentation results of the 'March' and the 'Jump' data are representatives of the local dynamic behaviors of different mesh regions. As can be seen in Figure 5, our segmentation model can not only automatically determine the number of segments, but also divide the mesh based on the local movements and group the disconnected regions with similar behaviors. • The 'Cloth' animation in Figure 5(II) is firstly segmented into 4 different highly deformed regions while dropping onto the table. Then, our approach generated 3 segments, i.e., 2 waving corner regions with deformed wrinkles and 1 relatively static region. • From the segmentation results of the 'Horse' animation in Fig-

ure 5(III), we can observe the 4 leдs are classified into the same group when moving towards the same direction; otherwise, they form different spatial groups. Similarly, the 'tail' is grouped with the 'trunk' region in the case of the absence of distinct movements, or it is divided into two groups if bended.

Parallel computing. As presented in Section 4.3, the actual temporal segmentation is applied to each spatial segment independently, which can be accelerated through parallel computing. In this experiment, we implement the actual temporal segmentation step with parallel computing on 4 cells. The computational time is shown in the column 's p ' in Table 1. Compared to the single thread implementation (column 's' in Table 1), the average efficiency has been improved by 10.71%, while it can be improved even up to 35.48% ('Cloth' data). It is noteworthy that the decompression time is less than 0.3s for 'March', 'Jump' and 'Handstand', less than 0.06s for 'Horse' and 'Cloth', and less than 0.65s for 'Flag'. This is important for applications that require a fast decompression such as bandwidth-limited animation rendering and display.

Compression results. Table 1 shows the different configurations of our spatio-temporal segmentation approach for the compression of the experimental data (µ = 0.99). For each of the data with different parameters, we highlight the best 'Rate', 'KGError', and 'Timing' in bold fonts. We present and discuss the compression results in reference to the following parameters:

• ϵ. It is a smoothing parameter for the initial temporal segmentation (Section 4.1). This parameter can be empirically chosen based on the target frame rate and the mesh complexity. • γ init . If we increase γ init for the initial temporal segmentation, the computing time may be significantly increased since the time complexity of the initial temporal segmentation (Section 4.1) is O(|γ init | 2 ). On the other hand, its influence on KGError is limited.

Moreover, bpv f tends to decrease for most of the experimental data (except the 'Handstand' data). • γ act . As can be seen in Table 1, the change of γ act does not significantly affect any of bpv f , KGError, and the computing time. This is because most of the actual temporal segmentation boundaries are found before reaching γ act . • ω. By increasing ω from 4 to 8, we do not observe the significant changes of the evaluation metrics. This is because our 2-stages spatial segmentation can automatically converge to the optimal number of spatial segments. Moreover, the multi-thread implementation of our approach significantly improves the computational efficiency (see the 'Timing' column in Table 1). Therefore, in general ω tends to be set to a small number. In fact, based on the previous studies [START_REF] Karni | Compression of soft-body animation sequences[END_REF][START_REF] Luo | Compression of 3D mesh sequences by temporal segmentation[END_REF], ω cannot be a big number because the bit rate will increase sharply due to the additional groups' basis. In our experiments, we empirically set ω = 4 because our experimental computer has a CPU of 4 cells.

Comparative Studies

We also compared our method with the method in [START_REF] Sattler | Simple and efficient compression of animation sequences[END_REF]] (called as and the 'Original Simple' method in this writing), which is a non-sequential processing compression method. Additionally, we adopted the idea in [START_REF] Aris | Adaptive compression of animated meshes by exploiting orthogonal iterations[END_REF] which cuts an animation into temporal blocks of the same size. Then, we can simulate the sequential processing of the existing compression methods, including the 'Original Soft' in [START_REF] Karni | Compression of soft-body animation sequences[END_REF]] and the PCA-based methods, to compress each block in order. We call the adapted approaches as the 'Adapted Soft' and the 'Adapted PCA'.

In order to make fair comparisons, the block size of the adapted methods is approximately set to the average of |ατ | for each of the experimental data. Note that we have not included an 'Adapted Simple' method, which can be obtained by similarly adapting the 'Original Simple' method, into the comparison due to the extremely high computational cost of the 'Original Simple' method [START_REF] Sattler | Simple and efficient compression of animation sequences[END_REF], which is unsuitable for sequential processing.

KG Error versus bpvf. Figure 6 shows the comparisons of an example between our method and the other methods. As can be seen from this figure, our method shows a significantly better performance than the adapted methods. That is, with the same bpv f in the range of [2, 6.5], our method can always reconstruct the 'Cloth' animation with a much smaller KG Error. Note that the 'Original Soft' method has a better performance when bpv f < 2. This is because the 'Cloth' data contains a large portion of nearly static poses, which means the animation has significant temporal redundancies. Thus, the non-sequential precessing method ('Original Soft') takes this advantage by treating the entire animation. However, our method runs much more efficiently: on average, 14.5 seconds consumed by our method, 32.5 seconds consumed by the 'Original Soft' method, and 4421.9 seconds consumed by the 'Original Simple' method. Moreover, our method also provides a fine option for users who prefer high qualities after compression with slightly more storage cost, e.g., bpv f > 2.

Reconstruction errors. Figure 7 shows the heat-map visualizations of the reconstruction errors by our method and the other methods. Overall, our method can achieve smaller reconstruction errors with lower bpvfs for the experimental data. We describe the comparative results in details as follows:

• Comparisons with the adapted methods. As can be seen in the left and the middle of Figure 7, high reconstruction errors occur randomly on the mesh using the 'Adapted Soft' method, as it is based on the linear prediction coding, which does not explicitly constrain the spatial affinities. For the 'Adapted PCA' method, high reconstruction errors occur in the regions of the vertices with rapid movements. • Comparisons with the non-sequential processing methods. As can be seen in the right of Figure 7, for the 'Cloth' data, the 'Original Soft' method behaves with similar symptoms as the 'Adapted Soft' method. The 'Original Simple' method returns high reconstruction errors on the table-top region because this method groups the vertex trajectories based on the entire mesh sequence, which constrains neither the temporal affinities in the local temporal subsequences nor the local spatial affinities. In addition, our method is significantly faster than the 'Original Soft' method and the 'Original Simple' method: our method consumed 17.78 seconds, the 'Original Soft' consumed 36.88 seconds, and the 'Original Simple' consumed 4622.36 seconds. • Our method. Based on the above findings, our method avoids local extreme reconstruction errors using the specially-designed spatio-temporal segmentation to exploit both the spatial and the temporal redundancies. This advantage becomes more significant when periodically dynamic behaviors either spatially or temporally occur in the animation. In addition, our method runs much more efficiently, compared to the non-sequential methods (i.e., 'Original Soft' and 'Original Simple').

Limitations

The main limitation of our current model is the configuration of the parameters needed for the spatio-temporal segmentation scheme.

To investigate this issue, we have conducted experimental analysis on the parameters in Section 5.2. Based on our analysis, the tuning of the parameters only has limited influence on the compression results. Using the 'Horse' data in Table 1 as an example, the compression does not change when we modify γ act from 20 to 30. This is because our approach often detects a temporal segmentation boundary before reaching γ act , case (II) in Section 4.5.

Another limitation of our model is the computational cost. Although we have implemented some parts of our spatio-temporal segmentation model through parallel computing and its computational time is superior to those of the existing non-sequential processing based compression methods, it requires further design for a frame-by-frame segmentation update scheme towards the real-time compression of 3D mesh animations in the future.

Figure 1 :

 1 Figure 1: Pipeline overview of our spatio-temporal segmentation scheme for compression. (I, II, III, and IV) are the 4 types of the segmented animation blocks, which are explained in Section 4.5. Note that b denotes the length of an initial/actual temporal segmentation, γ init and γ act are the maximal possible lengths for the Initial Temporal Cut and the Temporal Segmentation (short for the Actual Temporal Segmentation), respectively.

Figure 2 :

 2 Figure 2: An example of the initial temporal segmentation of the 'March' data, with the pairwise frame based kCCA matrix (Eq. 3) in the top panel and the MMN curve (Eq. 2) in the bottom panel. b is the detected boundary frame. The colorbar indicates the small (blue) and large (red) kernels.

  ; Wuhrer and Brunton 2010], we compute the Maximal Edge-length Change (MEC) within |τ | frames for each vertex pair, see Figure 3(I). (2) Binary labeling of vertices. We fit the MEC of all the edges as an exponential distribution epd, see the top of Figure 3(I).

Figure 3 :

 3 Figure 3: Pipeline of the vertex clustering within an initial temporal cut of the 'March' data: (I) Maximal Edge-length Change (MEC) for all the edge pairs and their distributions, (II) binary labeling of vertices, (III) the rigid clusters resulted from the initial vertex clustering, and (IV) the rigid cluster grouping results.

Figure 4 :

 4 Figure 4: Illustration of the temporal segmentation. The top row shows a sampled mesh sequence, with a bounding box as a sliding window. The size of the window is dynamically determined as the length of the initial temporal cut, i.e., |τ |. The bottom-left shows the vertex grouping of the initial temporal segment, and the bottom-right contains the change of the number of PCs for each vertex group in the sliding window. γ act is the maximal possible delay, and |ατ | is the detected temporal segmentation boundary.

Figure 5 :

 5 Figure 5: The spatio-temporal segmentation results of the experimental data: (I)'March', (II)'Cloth', (III)'Horse', (IV)'Jump'. Note that colors only indicate the intra-segment (not inter-segment) disparities. See more results in the supplemental materials.

Figure 6 :

 6 Figure 6: Comparisons on the 'Cloth' animation between our model (ω = 4, γ init = 40, γ act = 100) and the 'Adapted Soft' (block size = 100), the 'Adapted PCA' (block size = 100), and the 'Original Soft'.

ACKNOWLEDGMENTS

This work has been in part supported by the National Natural Science Foundation of China (No.61602222, 61732015, 61762050, 61602221), the Natural Science Foundation of Jiangxi Province (No.20171BAB212011) and the Key Research and the Development Program of Zhejiang Province (No. 2018C01090), and US NSF IIS-1524782.

Figure 7: The reconstruction errors of the compression by using our method, 'Adapted Soft', 'Adapted PCA', 'Original Soft' [START_REF] Karni | Compression of soft-body animation sequences[END_REF]] and the 'Original Simple' [START_REF] Sattler | Simple and efficient compression of animation sequences[END_REF]]. The colorbar indicates the reconstruction errors from low (blue) to high (red).

CONCLUSION

In this paper, we present a new 3D mesh animation compression model based on spatio-temporal segmentations. Our segmentation scheme utilizes a two-stages temporal segmentation and a twostages vertex clustering, which are greedy processes to exploit the temporal and spatial redundancies, respectively. The main advantage of our scheme is the automatic determination of the optimal number of temporal segments and the optimal number of vertex groups based on global motions and the local movements of input 3D mesh animations. That is, our segmentation scheme can automatically optimize the temporal redundancies and the spatial redundancies for compression. Our experiments on various animations demonstrated the effectiveness of our compression scheme. In the future, we would like to extend our spatio-temporal segmentation scheme to handle various motion representations, which can be potentially used for various motion-based animation searching, motion editing, and so on.