domain: info.info-ts info math

text:
Mariia Fedotenkova
email: mariia.fedotenkova@inria.fr

Axel Hutt

Research report: Comparison of different time-frequency representations

Chapter 1 Introduction

Time-frequency representation (TFR) is analysis methods which expands stationary spectral analysis to non-stationary cases. It is used for distinguishing different frequency components presented in the signal as well as their evolution in time.

This report compares different time-frequency representations based on previous studies (see [START_REF] El | A comparative study of time-frequency representations for fault detection in wind turbine[END_REF], [START_REF] Jones | A resolution comparison of several time-frequency representations[END_REF] and [START_REF] Thayaparan | Linear and quadratic time-frequency representations[END_REF]) and implementations (see [START_REF] Auger | Time-frequency toolbox. for use with MATLAB[END_REF]) but in application to EEG data analysis. This work compares several time-frequency representations taking into account speed, accuracy and localization of each method and chooses the best one for the particular problem. The methods discussed in this report are:

• Short-time Fourier transform and spectrogram (STFT);

• Continuous Wavelet transform and scalogram (CWT);

• Wigner-Ville distribution;

• Choi-Williams distribution;

• Reassigned spectrogram.

The comparison was held on two test data sets and afterwards the winning TFR was applied to the real data.

Chapter 2

Overview of time-frequency representations

First step of the analysis was held on the test data presented by sum of sines for all studied timefrequency representations.

Test data description

Test signal is a sum of oscillations with three different frequencies appearing at different times. It is presented on Fig. 2.1. This chapter describes the data used for the first step of the analysis and also gives an overview of methods used to obtain time-frequency representations. Original data consist of sum of three signals, which are solution of Lotka-Voltera equations with added high-amplitude oscillations on peak values and small noise oscillations elsewhere (σ = 0.1). It is of the form:

     z 1 = I 1 sin(ω 2 1 t) + ξ, ω 1 = 2π • 0.4 z 2 = I 2 sin(ω 2 2 t) + ξ, ω 2 = 2π • 0.6 z 3 = I 3 sin(ω 2 3 t) + ξ, ω 3 = 2π
• 1 where I is the term, which amplifies oscillations at peak values of the solution of Lotka-Voltera system; ξ is a noise term (σ = 0.1); ω is angular frequency (in rad). Frequency in Hz can be derived as follows:

sin(2πf

1 t) = sin(ω 2 1 t) = sin(2π(2π • 0.4 2 )t) ⇒ f 1 = 2π • 0.4 2 and      f 1 = 2π • 0.4 2 = 1.0053 f 2 = 2π • 0.6 2 = 2.2619 f 3 = 2π • 1 2 = 6.2832
Those frequencies can be seen from PSD of each signal on the fig. 2.2. With known frequency and times at which they appear the ideal time-frequency representation can be constructed (see fig. 2.3). In addition small amount of blur was added.

(a) z 1 (b) z 2 (c) z 3

Continuous wavelet transform and scalograms

Similarly with STFT and spectrogram, scalogram is squared value of continuous wavelet transform (CWT) ([1], [START_REF] Torrence | A practical guide to wavelet analysis[END_REF] ). On fig. 2.5, the left picture is a real Morlet wavelet, the middle one is 4 th derivative of complex Gaussian and the right one is a complex Morlet with center frequency 0.95 and bandwidth parameter 2. Complex wavelets give more smoothed results because they include both real and imaginary parts, which gives broader peaks. Pseudo WVD (smoothed in frequencies) reduces those terms, but still not enough (middle image). Smoothed pseudo WVD (smoothing of WVD both in time and frequency) gives good results but requires a lot of time for computations.

To obtain these representation analytical signal was used. It was obtained from original signal using hilbert MATLAB function. It was done due to the fact that with discrete WVD aliasing may appear at half of Nyquist frequency (in contrast to Fourier transform).

In addition, 1%-thresholding was applied to these TFRs. It was done, because WVD gives negative values (see fig. 2.9). In order to perform comparison with known TFR those negative components were truncated (see fig. This chapter provided overview of existing time-frequency methods and further chapter will summarize the quantitative analysis of the latter.

Chapter 3

Comparison of time-frequency representations

Quantitative comparison of all methods

To quantitatively compare different TFRs different approaches were used. The first method is to measure classification error as difference between estimated TFR and ideal representation. To do so, firstly, they were 0-1 normalized, where white color Next criterion is a mean squared error (MSE) [START_REF] Zhou | Mean squared error: Love it or leave it? a new look at signal fidelity measures[END_REF] between two TFRs, given in percents:

M SE = I es -I id 2 N * 100,
where I es is an estimated TFR, I id is the ideal TFR and N is the number of elements in TFR.

Three next parameters are taken from [START_REF] Rankine | A quantitative comparison of non-parametric time-frequency representations[END_REF]. They are:

• The two-dimensional correlation between images ρ;

• The Instantaneous frequency correlation IF ;

• The Time-frequency resolution measure res.

For the first two parameters 1.00 means best quality of estimated representation, while for the third one it is 0.00.

All of this criteria for each method are presented in Table 3.1. For scalogram FFT-based method was used, since it slightly faster than one based on convolution, also to perform the comparison 3%-threshold was used.

Comparison for different noise levels

The previous data had the standard deviation of noise σ = 0.1. The following data has higher noise level, but the same TFR.

Signal with frequency modulations

This section introduces new test data set and also uses for the analysis only two leading methods:

Scalogram and Reassigned spectrogram.

Test data

The first half of the analyzed signal consists of sum of two sine waves with equal amplitudes but different frequencies Chapter 4

Results for different noise levels

Real data analysis

Data description

Data is obtained from 6 different patients and recorded during routine anesthesia, before and after surgical incision. Each data file contains:

• 30-second recording 2.5 -2 min before incision;

• 30-second recording 2 -2.5 min after incision;

• 10-minute recording around incision (1-47Hz filtered). High peak at the beginning of incision because of the tap on the electrode.

Signals are obtained from prefrontal montage and sampled at 128 Hz.

Patients can be divided in two groups by:

1. anesthetic drug they received: propofol or desflurane;

2. whether they had dreams or not during operation.

Signals in time domain are presented in the following figure. Time is given in seconds, top axes contains the signal centered around incision, the bottom left axes presents 30 seconds before incision and the bottom right one -30 seconds after.

Results for real data

The time-frequency representations were obtained from the signals using method of spectrogram reassignment, it was calculated only for pre-and post-incision short intervals.

The analyzed frequencies were chosen to be from 0 Hz to 64 Hz (half of the sampling frequency) with a step of 0.1 Hz. But the signal was 1-47 Hz filtered and looking at the fig. 4.4 it can be seen that no activity is present over 30 Hz, because of that it would be more appropriate to use smaller frequency range (1-30 Hz) but with higher resolution (smaller step). Frequency step should be chosen not less than 0.05 Hz, because duration of the signal T should be more than 1/∆f . Duration of the signals is 30 s, that means that ∆f > 1/T = 0.033. The time-frequency representations obtained with reassigned spectrogram can be seen on the figures below. Next, time-frequency representation was averaged over five frequency bands:

• delta: 1 -4 Hz;

• theta: 4 -8 Hz;

• alpha: 8 -12 Hz;

• beta: 12 -20 Hz;

• gamma: > 20 Hz. Chapter 5

Conclusion

Wigner-Ville and Choi-Williams distributions introduce a lot of distortions due to the interference terms and negative components. They also requires a lot of computational time.

Spectrogram better preserve energy distribution (comparing to wavelet transform), but it has issues with time/frequency resolution and gives rather rough estimates.

Wavelet transform and scalograms based on it provide good results, but add artifacts on frequencies close to zero (convolution-based method) or close to the edges (FFT-based method). Also it distorts energy distribution: high frequencies have less power comparing to lower ones. Also the TFR is blurred and it gets stronger with rise of frequency.

Frequency reassignment appears to provide indeed precise results from visual inspection. Quantitative characteristics (N c , N f , N m ) might have low values because most of the power on the TFR looks "reduced" comparing to certain small regions. So these regions get value 1 after normalization and the rest of frequency components are ≈ 0.5, which can be a reason of the error, since ideal TFR assumes 1 for all of the components.

Final notes:

1. Spectrogram computed by Matlab might overcome Spectrogram computed with TFTB in certain criteria, but it uses less time points, which is bad for recurrence plots.

2. Spectrograms shows better results than WV or CW distributions for this test data, because frequency components are linear, and WVD are considered to be better for nonlinear frequency laws.

3. Spectrogram reassignment was chosen as a method in future work due to:

• good localization both in time and frequency;

• high noise resistance;

• least power distortions.

4. Methods which could also be studied:

• filtered spectrogram;

• synchrosqueezing;

• multitapers.

Figure 2 . 1 :

 21 Figure 2.1: Test signal final.

Figure 2 . 2 :

 22 Figure 2.2: Power spectrums of signals z 1 , z 2 , z 3 .

Figure 2 . 3 :

 23 Figure 2.3: Ideal TFR.

2. 2 Time-frequency distributions 1 .

 21 Short-time Fourier transform and spectrogram Short-time Fourier transform (STFT) is windowed fourier transform obtained with window sliding along the signal, thus providing a matrix (time × frequency) instead of a vector (frequency) obtained with a normal Fourier transform. Spectrogram is squared values of STFT. On fig.2.4, first two plots are obtained with Matlab routine. Left picture represents a power spectral density (PSD) in a logarithmic scale: 10 • log 10 (|P SD|).Plotting only absolute values of PSD provides better readability (middle and right figures). Right picture is obtained using Time-Frequency Toolbox (TFTB)[START_REF] Auger | Time-frequency toolbox. for use with MATLAB[END_REF] and gives more smoothed results, tough slightly more blurred. Spectrogram obtained with TFTB of the size 300 × 3000, while for TFR computed with MATLAB it is only 300 × 22. It happens because MATLAB returns vector of time as points where the spectrogram was computed (window location). In both cases Hanning window of length 256 was used.

Figure 2 . 4 :

 24 Figure 2.4: Spectrograms.

Figure 2 . 5 :Figure 2 . 6 :

 2526 Figure 2.5: Continuous Wavelet transform.

Figure 2 . 7 :

 27 Figure 2.7: Scalograms obtained with different techniques.

  2.8).

Figure 2 . 8 :

 28 Figure 2.8: Wigner-Ville distributions after eliminating negative components.

Figure 2 .

 2 Figure 2.9: Wigner-Ville distribution with negative components.

Figure 2 .

 2 Figure 2.10: Choi-Williams distribution.

Figure 3 . 1 :

 31 Figure 3.1: Spectrograms obtained with TFTB and its differences with the ideal TFR.

Figures 3 .

 3 2 and 3.7 shows the signals under study in a time space (σ = 0.3 and σ = 0.5 accordingly), figures 3.3 -3.6 shows estimated TFRs for the signal with σ = 0.3 and figures 3.8 -3.11 for σ = 0.5, finally, Tables 3.3 and 3.4 give estimation of performance measures for both signals.

3. 2 . 1 σ = 0. 3 Figure 3 . 2 :

 21332 Figure 3.2: Test signal with higher noise level (σ = 0.3).

Figure 3 . 3 :

 33 Figure 3.3: Spectrograms (σ = 0.3).

Figure 3 . 4 :

 34 Figure 3.4: FFT-based scalogram using complex Morlet wavelet (σ = 0.3).

Figure 3 . 5 :Figure 3 . 6 :

 3536 Figure 3.5: Wigner-Ville distributions (σ = 0.3).

Figure 3 . 7 :

 37 Figure 3.7: Test signal with higher noise level (σ = 0.5).

Figure 3 . 8 :

 38 Figure 3.8: Spectrograms (σ = 0.5).

Figure 3 .

 3 Figure 3.9: FFT-based scalogram using complex Morlet wavelet (σ = 0.5).

Figure 3 .

 3 Figure 3.10: Wigner-Ville distributions (σ = 0.5).

  (a) Choi-Williams distribution (b) Spectrogram reassignment

Figure 3 .

 3 Figure 3.11: Choi-Williams distribution and Spectrogram reassignment (σ = 0.5).

(3 and 5

 5 Hz) and the second part consists of sinusoidal frequency modulation, where frequency changes from 3 to 8 Hz in accordance with sine function. Then white noise with different standard deviations (σ = 0.1, 0.3, 0.5) was added to the signal forming three different test signals. This section, as mentioned above, considers only scalogram (obtained using FFT, complex Morlet) and reassigned spectrogram. The signals in the time domain with different noise levels and its corresponding TFR can be seen on fig.3.13 and fig.3.12. Reassigned spectrograms and scalograms for the signals with different noise levels can be seen on the fig.3.14 -3.16 and quantitative comparison in the tab.3.5.

Figure 3 .

 3 Figure 3.12: Ideal time-frequency representation for the test signal.

Figure 3 .

 3 Figure 3.13: The signals in the time domain with different noise levels.

Figure 3 .

 3 Figure 3.14: Time frequency representations for the signal with noise level σ = 0.1.

Figure 3 .

 3 Figure 3.15: Time frequency representations for the signal with noise level σ = 0.3.

  (a) Patient #1201. Propofol. Dreaming. (b) Patient #1203. Desflurane. No dreaming.

Figure 4 . 1 :

 41 Figure 4.1: Real data in the time domain. Intervals centered at, pre-and post-incision.

  (a) Patient #4002. Propofol. Dreaming. (b) Patient #4010. Desflurane. No dreaming.

Figure 4 . 2 :

 42 Figure 4.2: Real data in the time domain. Intervals centered at, pre-and post-incision (Continued).

  (a) Patient #4011. Propofol. No dreaming (b) Patient #4012. Desflurane. No dreaming.

Figure 4 . 3 :

 43 Figure 4.3: Real data in the time domain. Intervals centered at, pre-and post-incision (Continued).

  (a) Patient #1201. (b) Patient #1203.

Figure 4 . 4 :

 44 Figure 4.4: Reassigned spectrograms for pre-and post-incision data.

Figure 4 . 5 :

 45 Figure 4.5: Reassigned spectrograms for pre-and post-incision data (Continued).

Figure 4 . 6 :

 46 Figure 4.6: Reassigned spectrograms for pre-and post-incision data (Continued).

  The results are shown on fig. 4.7 (a) Patient #1201. (b) Patient #1203.

Figure 4 . 7 :

 47 Figure 4.7: Reassigned spectrograms for pre-and post-incision data. Average over frequency bands.

Figure 4 . 8 :

 48 Figure 4.8: Reassigned spectrograms for pre-and post-incision data. Average over frequency bands (Continued).

  (a) Patient #4011. (b) Patient #4012.

Figure 4 . 9 :

 49 Figure 4.9: Reassigned spectrograms for pre-and post-incision data. Average over frequency bands (Continued).

Table 3 .

 3 Also average computational time was measured for all of the representation (see Table3.2).

		1: Performance measures for each method
			N c	N f	N m M SE	ρ	IF	res
	Spectrogram (Matlab)	97.59 1.51 0.89 0.103 0.865 0.744 0.0002
	Spectrogram (TFTB)	97.21 1.63 1.16 0.096 0.865 0.715 0.0002
	Scalogram (C.Morlet, FFT) 91.21 7.70 1.08 0.360 0.605 0.360 0.0016
	WVD	77.8 20.56 1.64 0.215 0.662 0.212 0.0008
	pWVD	79.34 19.27 1.39 0.123 0.824 0.257 0.0009
	spWVD	97.01 1.63 1.36 0.088 0.888 0.716 0.0001
	CWD	96.30 2.33 1.36 0.089 0.885 0.629 0.0001
	Reass. Spectrogram	57.81 40.44 1.74 0.250 0.648 0.133 0.0001
		Table 3.2: Computational times for each method
	Spectrogram (Matlab) Spectrogram (TFTB) CWT (Morl.) CWT (C. Gaus.)	CWT (C. Morl.)
	0.160	0.103			0.228	0. 222	0.322
	Scalogram (Morl.)	Scalogram (C. Gaus.)		Scalogram (C. Morl.)	Scalogram (TFTB)
	0.233	0.229			0.329		6.007
	WVD	pWVD		spWVD	CWD	Reas. Spectro.
	0.118	0.119			16.280	16.282	6.849

Table 3 . 3

 33 82.27 16.68 1.05 0.104 0.857 0.309 0.0010 Scalogram (C.Morlet, FFT) 88.82 9.96 1.22 0.312 0.580 0.321 0.0013 WVD 58.98 39.42 1.6 0.348 0.496 0.149 0.0038 pWVD 57.93 40.69 1.37 0.276 0.631 0.176 0.0044

		N c	N f	N m M SE	ρ	IF	res
	Spectrogram (Matlab)	83.68 15.48 0.83 0.107 0.859 0.286 0.0013
	Spectrogram (TFTB)
	spWVD	87.57 11.07 1.34 0.093 0.876 0.347 0.0005
	CWD	78.37 20.26 1.36 0.098 0.867 0.258 0.0005
	Reass. Spectrogram	57.83 40.43 1.74 0.253 0.649 0.122 0.0002
	3.2.2 σ = 0.5

: Performance measures for each method with higher noise level (σ = 0.3)

Table 3 .

 3

	4: Performance measures for each method with higher noise level (σ = 0.5)
		N c	N f	N m M SE	ρ	IF	res
	Spectrogram (Matlab)	59.65 39.38 0.97 0.136 0.820 0.182 0.0039
	Spectrogram (TFTB)	57.58 41.19 1.23 0.121 0.827 0.197 0.0033
	Scalogram (C.Morlet, FFT) 73.57 25.15 1.27 0.327 0.549 0.212 0.0031
	WVD	53.26 45.15 1.59 0.679 0.320 0.138 0.0082
	pWVD	52.31 46.23 1.45 0.577 0.414 0.158 0.0091
	spWVD	67.68 30.89 1.43 0.121 0.845 0.208 0.0019
	CWD	62.39 36.15 1.46 0.134 0.817 0.190 0.0017
	Reass. Spectrogram	57.42 40.85 1.73 0.273 0.606 0.118 0.0003

Table 3 .

 3 As it can be seen from the figures and the table spectrogram reassignment surpasses scalogram in noise resistance, localization and introduces least distortions. Next chapter will show results of spectrogram reassignment on the real data.

		5: Performance measures for each method
		N c	N f	N m M SE	ρ	IF	res
	σ = 0.1	Reass. Spectrogram 59.37 40.13 0.50 0.268 0.797 0.103 0.0001 Scalogram 73.76 25.74 0.5 3.250 0.243 0.138 0.0133
	σ = 0.3	Reass. Spectrogram 58.69 40.81 0.50 0.286 0.769 0.100 0.0001 Scalogram 72.42 27.08 0.5 2.201 0.238 0.135 0.0107
	σ = 0.5	Reass. Spectrogram 58.44 41.06 0.50 0.352 0.740 0.096 0.0002 Scalogram 62.63 36.7 0.5 1.807 0.228 0.116 0.0116


Rodolphe Vaillant

Gaël Guennebaud

Loïc Barthe

Brian Wyvill

Marie-Paule Cani

Robust Iso-Surface Tracking for Interactive Character Skinning

Keywords: CR Categories: I, 3, 7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Animation; Skinning with contact, Iso-surface tracking

 [START_REF] Vaillant | Implicit skinning: Real-time skin deformation with contact modeling[END_REF]) which locally deforms the mesh according to some preset weights and (c) our new technique which automatically produces a plausible skin elasticity (notice how the belly button stretches). On the right, the Dana's knee is bent with an extreme rotation angle. (d) Implicit Skinning fails to handle deep self-intersections while (e) our technique allows large bending angles and the self-intersection at the knee is handled correctly.

Abstract

We present a novel approach to interactive character skinning, which is robust to extreme character movements, handles skin contacts and produces the effect of skin elasticity (sliding). Our approach builds on the idea of implicit skinning in which the character is approximated by a 3D scalar field and mesh-vertices are appropriately re-projected. Instead of being bound by an initial skinning solution used to initialize the shape at each time step, we use the skin mesh to directly track iso-surfaces of the field over time. Technical problems are two-fold: firstly, all contact surfaces generated between skin parts should be captured as iso-surfaces of the implicit field; secondly, the tracking method should capture elastic skin effects when the joints bend, and as the character returns to its rest shape, so the skin must follow. Our solutions include: new composition operators enabling blending effects and local self-contact between implicit surfaces, as well as a tangential relaxation scheme derived from the as-rigid-as possible energy to solve the tracking problem.

Introduction

Many computer graphics applications include the animation of synthetic 3D characters. One challenging aspect is the production of realistic skin deformations, especially at joints. To achieve realism, self-collisions must be avoided and contact surfaces generated instead. In addition, being able to add muscle inflation and skin stretching during body animation is highly desirable. Two main approaches have been used for skinning characters so far: geometric skinning and physically-based simulation. Geometric skinning techniques such as smooth skinning or dual quaternions enable realtime animation, but are not able to perform contact deformations nor reproduce plausible skin elasticity (i.e., its ability to stretch and slide over the character's volume). In contrast, physically-based simulation can be used to achieve these effects, however, long computation times prevent interactive feedback. In addition user control of the muscle setup (when necessary) can be very tedious.

Recently, skinning with contact deformations has been achieved in real time using the implicit skinning method [START_REF] Vaillant | Implicit skinning: Real-time skin deformation with contact modeling[END_REF]]: the character's volume is approximated by a set of local, 3D scalar fields (Figure 3 (b)-top). The scalar fields approximate the limbs, and at each animation frame, they are rigidly transformed by the skeleton and combined with an appropriate composition operator, producing a single scalar field (Figure 3 (b)-bottom). To obtain the final deformation, mesh vertices are deformed using standard geometric skinning and then iteratively projected onto their original iso-surface within the scalar field.

While the results of this technique are encouraging, its use of standard geometric skinning as an initial solution requires specifying skinning weights, which are tedious to define. Moreover, as observed by the authors, the success of the correction is conditioned by the quality of the initial skinning solution. When the initial vertices are too far from their ideal final position, which always happens for large bending angles, the following problems occur:

• Foldovers and triangle inversions are produced when some triangles deformed by geometric skinning are nearly par- • Bad projections occur when self-intersecting vertices cross the scalar field's medial axis. This projects vertices on the wrong side of the medial axis. (Figure 2(b)-bottom).

• Mesh distortions (e.g. excessive stretch) introduced by the geometric skinning step (Figure 2(a)) are not corrected during the projection on the implicit skin (2(b)).

• The projection is unstable due to the fact that contacts are detected as gradient discontinuities in the 3D scalar field, but not modeled explicitly as iso-surfaces of the field.

In practice, these flaws restrict the range of possible character poses, making the method less useful in a production environment.

Our goal is to design a production friendly skinning technique, which takes advantage of the desirable features of the implicit skinning framework, but avoids the pitfalls. The solution presented in this paper captures skin contact and elasticity, as well as muscle inflation and skin sliding effects. It also provides more automation and robustness than previous methods while still being applicable at interactive frame-rates.

More precisely, our contributions include:

• A new implicit skinning pipeline, free from the need of any initial skinning solution and therefore avoiding both tedious specification of skinning weights and dependence on the quality of a geometric solution (Section 3).

• New composition operators enabling us to approximate the shape of the character's skin during animation with an implicit skin that includes contact surfaces with a controllable depth, as depicted in Figure 3(e) (Section 4).

• A robust tracking method based on a linear relaxation energy (inspired from the as-rigid-as-possible energy) to update the skin mesh with a plausible tangential distribution of vertices, as shown in Figures 1(c) and 2(c) (Section 5).

• Easy extensions to advanced effects such as muscle inflation (Figure 10) and skin sliding, thanks to the clear decoupling between volumetric deformations handled by the implicit skin and stretching effects handled by the skin mesh (Section 6).

2 Related Work

Skinning methods

The most popular methods to animate the skin of a character are geometrical [START_REF] Magnenat-Thalmann | Modeling of bodies and clothes for virtual environments[END_REF]Kavan and Žára 2005;[START_REF] Kavan | Geometric skinning with approximate dual quaternion blending[END_REF][START_REF] Kavan | Elasticity-inspired deformers for character articulation[END_REF][START_REF] Kim | Bulging-free dual quaternion skinning[END_REF].

In particular, linear blend skinning and dual quaternion skinning are very fast to compute and relatively simple to implement. At run time, they directly weight the skeleton transformations before applying them onto the corresponding mesh parts. Although associating skinning weights to mesh vertices can be done automatically [START_REF] Lasa | Geodesic binding for degenerate character geometry using sparse voxelization[END_REF], the artist often has to manually tune them to improve the results. Moreover, these two methods tend to either increase or decrease the volume around joints, and interpenetrations are generated when different skin parts come in contact.

Shape interpolation schemes [START_REF] Lewis | Pose space deformation: a unified approach to shape interpolation and skeleton-driven deformation[END_REF]] and example-based methods [START_REF] Mohr | Building efficient, accurate character skins from examples[END_REF][START_REF] Wang | Real-time enveloping with rotational regression[END_REF][START_REF] Weber | Context-aware skeletal shape deformation[END_REF]] can also be used to skin a character in real-time. Unfortunately, capturing contact deformations is not very practical since it requires the inclusion of all possible contact situations in the example database. Moreover, skin shapes cannot be easily interpolated or extrapolated in the neighborhood of a contact, the deformation being highly non-linear.

A solution to get more realistic results is to use physically-based simulation, for instance based on anatomical knowledge [START_REF]Anatomically-based models for physical and geometric reconstruction of humans and other animals[END_REF][START_REF] Teran | Robust quasistatic finite elements and flesh simulation[END_REF][START_REF] Mcadams | Efficient elasticity for character skinning with contact and collisions[END_REF]]. The necessary anatomical data can be semi-automatically transferred between different characters [START_REF] Cani | Anatomy transfer[END_REF]], yet, the intensive computation required usually prevents interactive feedback.

In addition, user control of the muscle setup (when necessary) can be very tedious. Methods such as elasticity skinning [START_REF] Mcadams | Efficient elasticity for character skinning with contact and collisions[END_REF] enables the computation of skin squash and stretch directly on the mesh by solving underlying physical equations, however, several seconds per animation frame are required for visualizing the results. Alternatively, one can use a set of training poses [START_REF] Teng | Simulating articulated subspace self-contact[END_REF] to accelerate the simulation.

Recently, skinning with contact modeling has been achieved in realtime using the implicit skinning method [START_REF] Vaillant | Implicit skinning: Real-time skin deformation with contact modeling[END_REF]]. In addition to generating visually compelling shapes at joints, this method provides a variety of transitional effects, from smooth skin junctions at the elbow for low bending angles to local bulges when fingers articulate. These effects capture at low cost some of the dynamic behavior provided by physically-based approaches. Similarly to geometric skinning approaches, implicit skinning requires the specification of skinning weights. Indeed dual quaternion skinning is used at each frame as an initial solution. Another flaw is that the method does not provide any implicit representation (i.e. the 3D scalar fields) of contact surfaces between skin parts: skin vertices trying to reach their associated iso-surface are simply stopped when they cross gradient discontinuities of the scalar field, causing instabilities. Lastly, as already discussed, the correction fails for deep self-penetrations between skin parts.

In this work, we revisit implicit skinning to avoid the use of an initial geometric skinning solution. Our method relies on both a new composition operator enabling us to model contact surfaces within the implicit skin, and a new, robust implicit surface tracking method applied to mesh points. We therefore review the previous work in these two domains below.

Composition operators for contact modeling

Starting from 3D scalar fields reconstructing the different parts of the input skin mesh, implicit skinning can make use of a variety of field composition operators to create the implicit skin, i.e. the 0.5 iso-surface around which the mesh is repositioned: If a fold needs to be generated as soon as the joint bends, [START_REF] Ricci | Constructive Geometry for Computer Graphics[END_REF] max operator is applied. More advanced effects such as a bulge-in-contact, or a smooth blend relies on specific gradient-based composition operators [START_REF] Gourmel | A gradient-based implicit blend[END_REF]. Gradient-based operators output a scalar field value which depends on both the input fields and their local gradients. Although they enable capturing the desired outside shape, these different operators have a shortcoming: they do not model the contact surface between the two shape parts. While this does not create any visual artifact, contact surfaces being usually hidden, it prevents handling the projection of mesh vertices in contact regions in a robust manner, resulting in some vertices being pushed towards the wrong side of the shape, causing the instabilities we mentioned.

A contact operator is a composition operator able to define a contact surface (within the 0.5 iso-surface in our case) where the combined implicit primitives collide. Such operators were developed for long [START_REF] Cani | An implicit formulation for precise contact modeling between flexible solids[END_REF]] and were reformulated recently for more flexibility [de Groot et al. 2014]. With these operators, the the implicit primitives to be combined remain completely disjoint: they are separated by the 0.5 iso-surface modeling the contact surface in the collision region. Unfortunately, this solution is not feasible for blending two implicit primitives representing the same limb near a joint. We need to design composition operators able to merge the primitives into a single volume in the rest pose, while locally producing a contact surface whose extent should increase as the joint's bend increases (see Figure 3(e)). In this work, we have developed specific contact composition operators, enabling this advanced combination of implicit skin parts.

Implicit surface tracking

Tracking an iso-surface of a dynamic 3D scalar field with a mesh remains a challenging problem in the general case. Indeed, the lack of natural tangential parameterization of the field level-sets makes it difficult to maintain an even distribution of mesh vertices.

In the case of implicit skinning, each mesh vertex needs to track a slightly different iso-surface over the animated implicit skin, in order to maintain geometric details. [START_REF] Vaillant | Implicit skinning: Real-time skin deformation with contact modeling[END_REF] solved the tracking problem by starting at each frame of a geometric skinning deformation of the mesh, and then interleaved Newton iterations in the field gradient direction with tangential mesh relaxation steps based on mean-value coordinates [START_REF] Hormann | Mean value coordinates for arbitrary planar polygons[END_REF] to project each mesh vertex on its iso-surface of interest. While effective in most cases, this method is subject to fold-over and triangle inversions, projection failures and inappropriate mesh distortions, as already explained in Section 1.

Another approach is to perform tracking iteratively: then, instead of using some additional initial solution, the mesh restarts at each frame from its position at the previous step. Most mesh-based iterative tracking techniques perform an orthogonal projection of mesh vertices onto the iso-surface, combined with some tangential relaxation controlling the vertex distribution. The first solution was to consider mesh vertices as particles and apply a global optimization scheme in order to maintain their distribution [START_REF] Witkin | Using particles to sample and control implicit surfaces[END_REF]. Other relaxation schemes make use of mass-spring systems [START_REF] Rodrian | Dynamic triangulation of animated skeleton-based implicit surfaces[END_REF] or of other physically inspired energies [START_REF] Bouthors | Twinned meshes for dynamic triangulation of implicit surfaces[END_REF]. These solutions cannot easily be adapted to our problem, because the relaxation of a particle system would deform the initial skin mesh, which we would like to keep unchanged since it stores geometric details. Also, mass-spring systems are known to be subject to fold-over. Lastly, all these methods would require dynamic remeshing through edge split and collapse to be kept robust. This is to be avoided in our case, where we would like the skin mesh to come back close to its initial configuration when the character comes back to the rest pose.

Active contours [START_REF] Xu | Snakes, shapes, and gradient vector flow[END_REF][START_REF] Sonka | Image Processing, Analysis, and Machine Vision[END_REF]] are also well known iso-surface tracking techniques. They focus on the accurate approximation of a given iso-surface rather than trying to deform an existing mesh while guaranteeing the time coherence of the deformation. Tracking can also be performed by assuming a translational transformation to extract a tangential component of the unknown deformation velocity field [START_REF] Stam | On the velocity of an implicit surface[END_REF] using differential information. This approach does not prevent the use of a relaxation scheme, and it appears less attractive for skinning, which mostly involves rotational transformations.

Making iterative tracking possible in an implicit skinning framework requires the introduction of a technique able to deal with contact surfaces and minimize the distortions in the vertex distribution.

In particular, mesh distortions should be avoided if the character comes-back to rest pose after several transformations. This motivates the introduction of a new, dedicated tracking scheme.

Overview

Our skinning framework is decomposed in two distinct stages.

Presets: Starting with a mesh partitioned according to the bones of its animation skeleton (Figure 3 (a)), each mesh part i is approximated by the 0.5 iso-surface of a smooth 3D scalar field fi : R 3 → R (Figure 3 (b)-top). Scalar fields fi have compact support and vary in [0, 1] with fi(p) < 0.5 for points p outside the limb and fi(p) > 0.5 for points p inside. They are generated by fitting Hermite-Radial-Basis-Functions to a subset of the mesh vertices. A single scalar field f : R 3 → R, whose 0.5 iso-surface (denoted as implicit skin) approximates the shape of the skin, is then computed by combining the fields fi with specific composition operators (Figure 3 (b)-bottom). In order to preserve the details of the original mesh, each vertex v stores its own iso-value f (v) within the resulting scalar field at rest-pose, to be projected back on this iso-value during animation.

While most of these presets follow the lines of [START_REF] Vaillant | Implicit skinning: Real-time skin deformation with contact modeling[END_REF], they differ from the nature of the composition operator used for f : we define a new family of gradient-based composition operators, able to capture contact surfaces between different parts of the implicit skin and to generate a smooth field around it (Section 4). Note that our new composition is flexible, and therefore enables us to progressively increase the extent of the contact surface while a joint bends.

During animation: The scalar fields fi are transformed according to their associated bone, resulting in a transformed version of the composed field f that embeds the implicit skin (Figures 3(c,d,e)bottom). In these Figures, we show planar sections of the 3D scalar fields. Blue lines are sections of outer iso-surfaces, grey lines are the section of the implicit skin and red lines are sections of inner iso-surfaces. This color convention is used throughout the paper.

In contrast with original implicit skinning, our tracking technique uses the mesh position at the previous time step, to deform the skin mesh according to the composed field f : Each mesh vertex starts from its final position at the previous animation step (Figure 3(c)) and is projected on its saved iso-value f (v) in the deformed field f while maintaining plausible elastic tangential mesh deformations and generating contact (Figures 3(d,e)). This is robustly done thanks to a new tracking method that successfully combines rigid transformations, standard Newton iterations following the gradients of the field f , and dedicated tangential relaxations (Section 5). Finally, a pass of Laplacian smoothing around con- tact regions filters the sharp creases at the boundaries of contact surfaces.

(a) (b) (c) (d) (e)
The two next sections respectively detail our solutions to the implicit composition and to the mesh tracking problems.

Implicit skin with contact modeling

We need to define a single 3D scalar field f whose 0.5 iso-surface, the implicit skin, produces self-contact regions upon collision while merging implicit primitives around the joints in a smooth manner.

The generation of a smooth C 1 continuous scalar field is required for robust iso-surface tracking (Section 5) as the gradient directions around the implicit skin need to point towards the closest (or nearly closest) point on the implicit skin, to limit mesh distortion at the tracking step. We present new operators for composing f in Section 4.1, and we explain how they are efficiently computed while retaining smoothness in Section 4.2.

Composition operators with contact

We solve the problem of enabling skin volumes to blend while sharing local contact surfaces using two new binary operators, organized in a composition tree [START_REF] Wyvill | Extending the csg tree -warping, blending and boolean operations in an implicit surface modeling system[END_REF]]. The leaves of the tree are the field functions fi each transformed in the same way as its corresponding bone and the composition operators are the interior nodes. A binary composition operator is a function g : R 2 → R combining two scalar fields f1 and f2 in a new 3D scalar field f k (p) = g(f1(p), f2(p)). Thus, each node of the tree defines a 3D scalar field and the root is the scalar field f including the implicit skin. Also, in our application we use a special type of composition operators called gradient-based operators.

As mentioned by [START_REF] Vaillant | Implicit skinning: Real-time skin deformation with contact modeling[END_REF], gradient-based operators g d : R 2 → R are parameterized by a function d : R → R whose argument is the angle α between ∇f1 and ∇f2. The parameter d controls the interpolation between two types of compositions, allows the discrimination of regions that blend according to the gra-dient angle and the application of different types of compositions.

In the case of skinning, it is used to localize blending or bulge-incontact effects in the fold only.

In the case of skinning, the effect of these operators is simpler to understand by noticing that the angle α between gradients in the folding regions is often close to the joint rotation angle. In the rest pose, no rotation is applied and thus, no specific skinning effect is to be performed. In that case, we need the 0.5 iso-surface of the scalar fields fi to be combined with a union while the rest of the iso-surfaces are smoothly blended around the implicit skin. This type of operator is called a clean-union operator [START_REF] Pasko | Function representation in geometric modeling: concepts, implementation and applications[END_REF][START_REF] Bernhardt | Implicit blending revisited[END_REF]. When a joint such as the elbow bends between α = 0 and α = π/2, the skin does not crease and no contact iso-surface need be generated. The use of a clean-union operator produces a smooth field around the implicit skin. For larger bending angles, α > π/2 at an elbow, the contact iso-surface has to be created using our contact operator. In this case the larger the bending angle, the deeper the contact surface goes towards the joint. Finally, when the gradients point in opposite directions (α = π), it means two disjoint skin parts have come into contact, for instance an arm against the torso, and a full contact operator should be applied.

Our first operator g dc is parameterized by the controller dc for contact handling. It smoothly interpolates between a clean-union (no contact) when dc(α) = 0 and a full contact (contact of maximal length) when dc(α) = 1. Therefore g dc can control the depth of the contact. The behavior presented above leads us to the plot of dc(α) depicted in Figure 4(a).

Figures 5(a,b) illustrate the relation between (a) the plot of a contact composition operator g dc (f1, f2) and (b) its effect on the composition of two scalar fields f1 and f2 for an intermediate value of dc.

In Figure 5(a), the abscissa (resp. ordinates) are values of f1 (resp. f2) and vertical (resp. horizontal) lines represents its iso-surfaces.

The operator defines the way the iso-surfaces of f1 and f2 are combined to produce the resulting scalar field f k . Values of g dc are those of f k when f k (p) = g dc (f1(p), f2(p)). Here, g dc = 0.5 on the 0.5 iso-surface of f1 (in yellow) up to its intersection with the 0.5 iso-surface of f2 (in purple), and on the 0.5 iso-surface of f2 up to its intersection with the 0.5 iso-surface of f2. The result is the union of these two surfaces. In addition, g dc = 0.5 on the green line which starts at the intersection of the combined 0.5 isosurfaces and goes inside the fields f1 and f2 following f1 = f2. This part of g dc = 0.5 is the contact surface that is shown with the same color on f k in Figure 5 ). The bulge-in-contact is performed by the inflation of the 0.5 isolines of f1 and f2 in Figures 6(b,c). During the animation of a finger, the parameter d b (α) is set as illustrated in Figure 4(b). As we can see, in that case, the bulge starts at low bending angles. Results produced by this operator are illustrated in Figure 13.

In practice, the composition operator g dc is applied by default at each skeleton joint, but the user can select g d b if desired. The binary composition tree (where implicit primitives are the leaves and composition operators are the nodes) is then built as follows: We start by creating nodes that link pairs of primitives with bulge-incontact composition g d b . This forms the bottom of the tree. Then all the other nodes are added on top using gradient-based contact g dc , until all primitives are combined within a single tree. This specific composition enables us to achieve both the required bulge surfaces at fingers and the contacts we are looking for.

Bi-harmonic evaluation

The advantage of using scalar fields with compact support is that they are bounded in space and thus, can be stored in 3D grids enabling fast evaluation. It is the same for composition operators: we pre-compute and store them in 3D grids for the three entries (f1, f2, d).

Most advanced composition operators, such as those based on gradient blending, have very complex closed form equations, so some of them must be evaluated numerically, for instance, using binary search. Building new composition operators this way is a long and tedious task: one must first derive intricate equations for each new operator and find the appropriate numerical evaluation method. Fortunately, in the case of skinning, only C 1 continuity is required for composition operators. This enables us to propose a simpler and more general way of constructing C 1 operators g directly into a discrete grid.

Our solution is to use bi-harmonic interpolation of some constraints: that is, the function g must satisfy the fourth-order bi-Laplacian equation ∆ 2 g = 0, where ∆ is the Laplacian differential operator, while matching constraints at specific input values, as those listed in the previous section. This choice ensures a minimal continuity of order C 1 . We directly solve this partial differential equation for each 2D slice (f1, f2) of the grid using a standard finite difference discretization. We use Dirichlet constraints on the system boundaries. The grid boundaries constraints are set according to the specific boundary properties of composition operators expressed by [START_REF] Canezin | Adequate Inner Bound for Geometric Modeling with Compact Field Function[END_REF] (Equation 1) and inner grid constraints are placed along the profile curve defining the 0.5 iso-curve of g(f1, f2) (Equation 2). Formally, we set:

g(f1, 0) = f1 , g(0, f2) = f2 , g(f1, 1) = g(1, f2) = 1, (1)
and g(f1, f2) = 0.5 ∀(f1, f2) ∈ profile curve.

(2)

In addition to these Dirichlet boundary conditions, the normal derivatives have to be constrained to zero on the g(f1, 0) and g(0, f2) axes to ensure a C 1 continuity at the boundaries of g.

Even though bi-harmonic interpolation does not guarantee the monotonicity of the function variation between constraints, we found the practical results satisfactory in our case without having to resort to more involved quadratic programming solvers [START_REF] Jacobson | Bounded biharmonic weights for real-time deformation[END_REF][START_REF] Jacobson | Smooth shape-aware functions with controlled extrema[END_REF].

Our different operators are thus generated by constraining g = 0.5 along a profile curve. For the contact operator g dc , the profile curve is composed of three line segments (in purple, yellow and green in Figure 5). The purple and yellow lines defined as the linear interpolation of the points (0, 0.5), (0.5, 0.5) and (0.5, 0.5), (0.5, 0) are fixed in all slices (f1, f2) ∀d ∈ [0, 1]. The green line is the linear interpolation of the points (0.5, 0.5) and (0.5

+ d √ 2 , 0.5 + d √ 2
). Aside from the operator g dc , the profile curve is in general not aligned with grid vertices (e.g., bulge operator). We constrain the grid vertices of the intersected grid cells with the shortest signed distance to the profile curve.

The bulge-in-contact operator g d b is constructed following exactly the same procedure (see Figure 6). The only difference are the purple and yellow profile curves that are defined with cubic B-splines whose control points are interpolated with respect to d b , between an aligned position when d b = 0 and the maximal bulge when d b = 1.

Using this evaluation procedure, our composition operators are precomputed once for all in 3D grids and at run time, a simple texture fetch is required for each evaluation. More generally, this approach enables the design of free-form composition operators: one only needs to specify a profile curve to generate a new operator.

Iso-surface tracking

One of the main challenge of our general implicit-skinning approach is to ensure that the initial mesh properly follows the isosurface defined by the time-varying global scalar field f . Even though f is usually obtained through the composition of implicit primitives fi which are rigidly transformed, the underlying implicit skin exhibits severe non-linear deformations around the joints. User defined non-linear deformations might also be applied to the fi to mimic some complex effects (see section 6). Moreover, in our context the mesh embeds many spatially varying data such as relief details and texture layers. Therefore, the topology and density of the mesh cannot be arbitrarily changed, and each vertex must carefully track its position on the implicit surface to avoid unnatural distortions.

To this end, we propose to track the iso-surface incrementally from one animation step to the next using a three stage procedure:

1. Each vertex is moved by applying the transformation difference of its nearest bone, as defined by the initial mesh partitioning (Figure 3 (d)).

2. Vertices are projected onto their respective iso-value (derived from the rest pose) along the gradient direction using Newton iterations.

3. This yields a good initial guess on which we can apply a tangential relaxation scheme to account for the stretch of the implicit skin while reducing distortions. Since stages one and two are already described by [START_REF] Vaillant | Implicit skinning: Real-time skin deformation with contact modeling[END_REF], the rest of this section focuses on the last one.

Tangential relaxation energy

The relaxation scheme must reduce the distortions while maintaining the mesh on the underlying surface. On the one hand, incrementally applying the relaxation based on mean value coordinate from the implicit skinning method does not reduce distortions sufficiently and contact is not captured adequately as illustrated in Figure 12(b). On the other hand, doing an accurate physical simulation of the skin would be prohibitively expensive for our purpose. We thus derived a novel relaxation scheme inspired by Sorkine et al. [2007] as-rigid-as-possible (ARAP) energy:

E({pi}) = n i=1 j∈N (i) wij (pi -pj) -Ri(pi -pj) 2 , (3)
where pi are the vertex positions in rest pose, pi are the unknown final positions, Ri is a 3 × 3 rotation matrix, and N (i) denotes the one-ring neighborhood of vertex i. The weights wij are classically defined from the cotangent formula [START_REF] Meyer | Discrete differential-geometry operators for triangulated 2-manifolds[END_REF], i.e., wij = 1 2 (cot(αij) + cot(βij)), where αij, βij are the angles opposite to the edge (i, j).

As motivated and detailed below, in our relaxation scheme the rotations Ri are prescribed a priori taking advantage of the specificity of skinning deformations. Indeed, the energy E is by construction invariant to translation while the rotation Ri locally cancels rigid transformations. In most applications of this kind of energy, including iso-surface tracking [START_REF] Bouthors | Twinned meshes for dynamic triangulation of implicit surfaces[END_REF], Ri is computed from the polar decomposition of the Jacobian of the unknown deformation. In addition to being an expensive process, this strategy exhibits several drawbacks. Firstly, this inserts a non-linear component in the energy, E which thus requires additional iterations to be minimized while introducing local minima. In addition, as illustrated in Figure 7(a,b), there is no smoothness constraints on the target rotations Ri and the result is only C 0 at the constraint boundaries. Moreover, we observed during animation that vertices appear to slide over the implicit skin instead of behaving as if the skin had elasticity, resulting in an unrealistic animation of the skin. Finally, as depicted in Figure 7(d), this strategy is well known to be subject to triangle inversions which are very tedious to avoid, even on a planar domain [Schüller et al. 2013].

We address all these drawbacks by fixing a priori the target rotations Ri. In order to leverage a coherent behavior over time and to recover the rest-pose, these rotations must not depend on the positions of the mesh vertices of the current animation step but only on the rest pose and affine transformations of the bones.

Following the first stage of our tracking procedure, a first choice consists of using the rotation associated with the nearest bone of vertex i for Ri. As illustrated in the first row of Figure 8, this strategy yields promising results despite the discontinuous nature of the underlying rotation field. An effective solution thus consists of performing a blend of the bone transformations as in standard geometric skinning methods. In this research we use the dual-quaternion interpolation scheme [START_REF] Kavan | Geometric skinning with approximate dual quaternion blending[END_REF]] with default harmonic weights [Baran and Popović 2007]. As discussed in more details in Section 6 and illustrated in Figure 8, we emphasize that both the choice of the blending technique and weights have only subtle effects on the final distribution of the vertices, and they do not affect the final shape or robustness of the tracking. Figures 7(c,e) show the superiority of our approach that better localize continuous tangential deformations and prevents triangle inversions.

On-surface constraints

The energy E has to be minimized while maintaining the vertices to their own field value. This is a non-linear constraint which implies an iterative scheme in which pi is obtained through a sequence p0 i , p1 i , . . . , where p0 i corresponds to the initial projection onto its respective iso-value, as obtained through the aforementioned second stage of the tracking system. At each iteration of the minimization, we first locally linearize this constraint by enforcing each vertex to lie on its respective current tangent plane. This can be accomplished by writing the unknown position pk+1 i as:

pk+1 i = pk i + u k i v k i u k+1 i v k+1 i , (4)
where u k i , v k i are two arbitrary tangent vectors of the current isosurface at pk i , and u k+1 i , v k+1 i are the 2D coordinates of pk+1 i in this local frame. Minimizing E under this constraint is still a linear and very sparse problem. Since there is no need to completely solve this problem before projecting the solution on the surface and iterating, we interleave one Jacobi iteration minimizing E with one Newton iteration projecting pk+1 i onto its respective field value. At each iteration, an approximation of the local tangent frame is obtained for free from the gradient computed during the Newton iteration.

The choice for Jacobi iterations is motivated by its simplicity to be parallelized on the GPU. Higher performance may be achieved using the conjugate-gradient method which can also be efficiently implemented on GPUs [START_REF] Weber | Efficient gpu data structures and methods to solve sparse linear systems in dynamics applications[END_REF]. In order to guarantee the availability of a very good initial guess p0

i and to avoid projection issues, the skeleton animation between two successive frames is decomposed into several virtual frames such that the maximal rotation angle of a bone is smaller than a third of the current bending angle. The intermediate skeleton configurations are linearly interpolated. Such a decomposition has little impact on the performance because a smaller animation step leads to a faster convergence.

Implementation and results

Implementation

In our prototype implementation, the deformation part of our skinning system is entirely executed on the GPU using CUDA. Parallelization is performed on a per-vertex basis with one kernel responsible for the computation of the initial guesses (rigid transformations and Newton projections) and a second kernel performing one iteration of the minimization (one Jacobi iteration, plus one Newton projection). This second kernel is applied multiple times per animation step until convergence. To accelerate the evaluation of the scalar field f , the scalar fields fi are indexed in a 32 3 grid, and both our composition operators g k and scalar fields fi are stored in 3D textures and evaluated using hardware trilinear interpolation.

Results were produced on an Intel Core i5-3570K at 3.4GHz, with 8GB of memory and a Geforce 670 GTX. The memory consumption is exactly the same as for the implicit skinning method: for the models presented in this paper, from 10 to 30MB of extra texture memory are required to store the operators and field functions.

During an interactive session, joints are usually adjusted one at a time. In this case, deformations occur mostly locally: only a small fraction of vertices have to be projected and relaxed, hence the density of the mesh has a small impact on performance. For a single joint, we thus observed a rate between 30 to 300 frames per second. More detailed statistics and performance results are given in Table 1 for complete animations and different time steps. As expected, the number of minimization iterations and thus the computation time highly depends on the magnitude of the skeleton transformation between two animation steps.

Results

Provided a mesh with an animation skeleton, our system automatically produces a default solution: partitioning and field-functions fi are automatically computed as described by Vaillant et al. [START_REF] Vaillant | Implicit skinning: Real-time skin deformation with contact modeling[END_REF], and skin elasticity and self-contacts are automatically and globally resolved during the animation.

As demonstrated in Figure 9, our approach produces satisfactory behavior without requiring any tedious tuning of skinning weights, while previous work still exhibit distortions even after an artist edits Dual Quaternion Skinning Implicit Skinning Our system the weights manually. The sensitivity of our method with respect to skinning weights is further evaluated and compared on a simpler configuration in Figure 8. Recall that in our system such weights are only used to blend the target rotation in the minimization of our tangential energy. Therefore, even large changes in the weights leads only to small differences in the distribution of the vertices, and the shape of the implicit skin remains unchanged. In contrast, the projection method of the implicit skinning technique fails if the weights are not smooth enough, hence, further editing of the skinning weights is often needed. The global nature of the skin sliding effect produced by our tangential relaxation is especially visible on the character's belly of Figure 1-left. In contrast to geometric or implicit skinning, with our method the belly button is stretched as the torso twists. As shown in Figure 11-top, this skin elasticity can be easily controlled by painting regions where the tangential relaxation is not wanted. In such regions, vertices are only rigidly transformed according to their nearest bone.

Global contacts can be observed in Figures 11(d) and 14. The shape of the deformation can also be adjusted by choosing between different types of composition operator to enable, for instance, subtle bulge in contact where this is desired as on the finger of Figure 13. More advanced deformations can also be performed by designing and controlling additional field functions. For instance, in Figure 10 the user added an extra implicit primitive with its associated bone to mimic the contraction of a bicep by applying a scaling proportional to the arm bending angle. Such a manipulation would not be possible without a proper tracking system taking care of the skin elasticity and contacts.

In addition to all these novel possibilities, our system permits reaching extreme bending angle without introducing artifacts. This can be seen in Figure 1(e) for a bending angle of around 150 degrees. This high robustness is due to the combination of both a proper modeling of surface contacts and the incremental nature of our novel surface tracking mechanism. Figure 12(a) shows that with a standard union operator, implicit skinning exhibits artifacts if it is implemented with skinning weights that are not smooth enough. In Figure 12(b) our novel contact operator, used alone in implicit skinning (i.e., without our novel tracking system), improves the projection artifacts but does not reduce mesh distortions. Finally, our novel tracking algorithm (Figure 12(c)) enables even larger bending Our system with our Our system with Implicit skinning default contact operator bulge and contact with bulge angles, a proper tracking of the fold depth and a reduction of mesh distortions.

Limitations

Owing to the incremental nature of our tracking technique, our approach gets closer to a physical simulation of the skin than the original implicit skinning method. Even though our approach is considerably simpler and faster than physical simulations, it still has some of their drawbacks. Firstly, the result for a given animation step is not entirely determined by the skeleton pose. Because of the non-linearity introduced by the projections on the implicit skin the result might depend on the complete history of the animation. For a similar reason, the result might also depend on the chosen time step.

As with any simulation, choosing the right time step is a trade-off between the computation cost and the numerical stability: a very large time step will lead to very poor initial guesses with similar problems as implicit skinning, while a very small time step will be prohibitively expensive.

In practice, as long as no incorrect projection occurs during animation, our solution is almost insensitive to a particular time step. The automatic solution proposed in Section 5.2 to choose an appropriate time step strives to make sure that the initial guess obtained through rigid transformation of the previous animation-step does not lead to a projection in the wrong direction. The proposed solution is rather conservative as, for instance, we observed that bending an arm from a fully open to a π/2 angle in a single time-step is correctly handled by our tracking system. Higher performance could thus be achieved Implicit Skinning Our complete system

Conclusion

In this work, we have presented a new skinning technique that takes advantage of the best features of implicit skinning, and makes the method more robust and suitable for production pipelines. We showed how to properly model and control contact surfaces using advanced implicit composition operators. To this end, we designed a process for the creation of free-form composition operators based on biharmonic interpolation over a discrete grid. The scope of this technical contribution goes beyond the context of implicit skinning, as this process enables interactive design of custom composition operators which is a central topic in implicit modeling. We also showed how to track a deformable implicit surface by deriving a linearized relaxation energy from Sorkine et al.'s as-rigid-as-possible energy with some a priori knowledge on the expected deformation that we can extract from the animation skeleton. As a result, our skinning system is robust, handles contact even when surfaces are deeply interpenetrating, and maintains a plausible solution where implicit skinning fails. The method thus reduces the amount of necessary user interaction, while maintaining interactive performance, and enables new features such as skin elasticity.

Based on these well controlled skin deformations, several new avenues could be explored. For instance, very fine control of the tangential deformations could be studied by setting a few vertices in key-frames. Also, performances could be improved by investigating more efficient GPU solvers and better automatic adjustment of the time steps. Automatically computing the scalar fields at complex joints is still a challenging problem which deserves some future specific attention. Finally, the more general use of our technique on deformable dynamic 3D objects is still to be explored.

Figure 1 :

 1 Figure 1: (a) The Jeff model in rest pose. Its shoulders are rotated and skinned with (b) implicit skinning (by[START_REF] Vaillant | Implicit skinning: Real-time skin deformation with contact modeling[END_REF]) which locally deforms the mesh according to some preset weights and (c) our new technique which automatically produces a plausible skin elasticity (notice how the belly button stretches). On the right, the Dana's knee is bent with an extreme rotation angle. (d) Implicit Skinning fails to handle deep self-intersections while (e) our technique allows large bending angles and the self-intersection at the knee is handled correctly.

Figure 2 :

 2 Figure 2: Illustration of different skinning solutions on a cylinder deformed with a single joint. (a) Dual quaternions produce bulge and self-intersection artifacts. (b-top) The implicit skinning corrects these artifacts, (b-bottom) but the projection stage fails for too large bending angles. (c) Our skinning technique avoids selfintersections, produces a contact surface and minimizes mesh distortions, even for very large bending angles.

Figure 3 :

 3 Figure 3: Presets: (a) Initial mesh partitioned according to its skeleton. (b) (top) Each partition is approximated by a scalar field fi and (bottom) composed with our new operators to produce the implicit skin. (c) (d) (e) We illustrate one step of our animation algorithm (top row the mesh deformation, bottom row the isolines of the implicit skin f ). During animation, (c) the deformed mesh parts at an initial step n (d) are rigidly transformed at the step n + 1. Then, vertices are projected onto their own iso-value f (v) (red arrows) and tangentially relaxed (blue arrows). (e) The final result with the contact surface illustrated on the mesh (top) and in grey in the implicit skin (bottom).

Figure 4 :

 4 Figure 4: (a) Plot of dc(α) to get contact at joints, and (b) d b (α) to mimic the inflation of tissues in contact situation.

  (b). This line always starts at the intersection of the combined primitives and its length is controlled by dc. Other iso-lines of g dc smoothly link values of f1 to values of f2 by following the shape of g d = 0.5, resulting in a smooth distance field around g dc = 0.5 in Figure5(a) and thus around f k = 0.5 in Figure5(b).By designing the same operator, but without contact (i.e., when dc = 0), we obtain the clean-union operator presented in Figure5(c) whose effect on scalar fields is depicted in Figure3(c). On the other hand, with a contact line of maximal length (i.e., when dc = 1), we get the full contact operator shown in Figure5(d).Our second operator g d b is built similarly. It interpolates both the depth of the contact and the inflation of the bulge-in-contact between a clean union operator when d b (α) = 0 (no contact and no bulge, Figure6(a)) and a full contact operator with maximal bulge when d b (α) = 1 (Figure 6(c)). Figure 6(b) illustrates this operator with intermediate contact length and bulge-in-contact (d b (α) = 1 2

Figure 5 :Figure 6 :

 56 Figure 5: (a) Our contact operator g dc (f1, f2) for an intermediate value of dc and (b) its application to the composition of two field functions articulated by a joint. In all figures, the contact surface is in green. (c) Clean-union operator (dc = 0, no contact) and (d) full contact operator (dc = 1).

Figure 7 :

 7 Figure 7: Standard ARAP energy (a), (b), (d) is compared to our dual-quaternion based rotation targets (c) and (e). In (a) only the extremities of the tube are fixed while in (b) a larger part is fixed. Both cases lead to C 0 deformations on the constraints with an unplausible distribution of vertices inside the domain. In (d) triangle inversions occurred. In (c) and (e) these problems are avoided and our technique produces an adequate result.

Figure 8 :

 8 Figure 8: Comparison of the sensitivity of the vertex distribution regarding variations of skinning weights (shown in color code) for the Implicit skinning (left) and our method (right). From top to bottom, skinning weights are rigid, smooth and smoother respectively.

Figure 9 :

 9 Figure 9: Comparison of different skinning method using manually edited skinning weights (top), and automatic weights (bottom).

Figure 10 :

 10 Figure 10: Reproduction of the inflation of a bicep by the addition of an extra primitive (a) and its animation using our system (b). Owing to the prominent shape of the muscles, implicit skinning yields projection artifacts (c).

Figure 11 :

 11 Figure 11: The Armadillo's knee (a) is bent using (b) our default solution, and (c) by disabling the tangential relaxation for red areas. (d) Global self-contact between the belly and the leg.

  Figure 12: (a) Without proper modeling of the surface contact, vertices tend to be projected towards the exterior of the contact region thus leading to incorrect projections. (b) Our novel operators address this precise issue, unfortunately, tracking with implicit skinning arbitrarily extends the contact region in an uncontrollable manner. (c) Our system adequately tracks the skin fold.

Figure 13 :

 13 Figure 13: Comparison of our bulge with contact operator.

Figure 14 :

 14 Figure14: Comparison of our system to implicit skinning on the model shown in rest pose in Figure1a. The arm collides too deeply with the belly to be properly recovered by implicit skinning while the chin is unrealistically creased. These issues are properly addressed by our system.

Table 1 :

 1 Statistics and performance of our system and implicit skinning for various models and time steps when all joints are animated simultaneously.

Acknowledgments

This work was partially funded by the IM&M project (ANR-11-JS02-007) and by the advanced grant EXPRESSIVE from the European Research council (ERC-2011-ADG-20110209). Partial funding also came from the Natural Sciences and Engineering Research Council of Canada, the GRAND NCE, Canada and Intel Corps. Finally, this work received partial support from the Royal Society Wolfson Research Merit Award. We wish to thank Laura Paiardini from Inria Grenoble for providing us with the very nice Armadillo's skeleton animation and Garrett Pond from Brigham Young University for the modeling of the Jeff model (the big guy).


Claire Lartigue
email: lartigue@lurpa.ens-cachan.fr

Yann Quinsat
email: quinsat@lurpa.ens-cachan.fr

Charyar Mehdi-Souzani

Zuquete Alexandre

Shadan Guarato

Yann Quinsat

Alexandre Zuquete-Guarato
email: alexandre.zuquete-guarato@renault.com

Shadan Tabibian
email: shadan.tabibian@renault.com

Voxel-based path planning for 3D scanning of mechanical parts

Keywords: Scanning path, visibility, Voxel-map, digitizing system 1

come

INTRODUCTION

Within the context of 3D surface scanning, path planning is still a challenge to obtain a complete representation of parts in a minimum time. More generally, path planning is defined as a set of ordered points of view allowing the part digitizing without collision. Each point of view corresponds to one sensor configuration relatively to the part and allows the digitizing of one portion of the part surface. As the common objective is to minimize cycle time, the main difficulty linked to scan path planning is to define a sensor trajectory, free from collision, that leads to a good surface representation, i.e. slightly noisy and complete. This issue can be expressed as the minimization of the number of points of view while respecting quality criteria, and avoiding collisions and occlusions. A point of view (,) pMd r is defined by a point M, and a direction d r

, and accounts for the sensor positioning (distance and orientation) relatively to the part surface (Fig. 1).

Finding the optimal path planning is an issue widely addressed in the literature, in particular for parts defined by a CAD model. Most of the methods are based on the concept of visibility. The visibility can be considered either from the part standpoint -what are the directions from which a point of the part surface is visible -or from the sensor standpoint -what is the surface area visible from the sensor point of view. Visibility is generally linked with the local normal at the point. Tarbox and Gottschlich [START_REF] Tarbox | Planning for complete sensor coverage in inspection[END_REF] propose to plan the path of a laser sensor by using the concept of measurability matrix. The measurability matrix results from the intersection of two matrices, one corresponding to the laser visibility the second one defining the CCD visibility. Starting from a CAD model, Xi and Shu determine optimal digitizing parameters so that the digitized surface is maximized. They split the surface into sections, each one defining a surface portion. For each section, the optimal positioning is obtained by aligning the field of view with the projected surface frontier [START_REF] Xi | CAD-based path planning for 3-D line laser scanning[END_REF]. Son et al. propose a method for path planning dedicated to a laser-plane sensor mounted on a CMM. The CAD model is first divided into functional surfaces. Each surface is thus sampled into points and corresponding normal vectors, with the objective of finding the minimum number of sensor configurations. The sampling is constrained by the size of the laser beam so that the distance between two consecutive points is less than the beam width. These two consecutive points can be digitized with the same sensor configuration if the angle between the normal vectors is less than twice the view angle of the sensor [START_REF] Son | Automated laser scanning system for reverse engineering and inspection[END_REF]. Within the context of 3D inspection Prieto et al. [START_REF] Prieto | Accuracy Improvement of contactless sensor for dimensional inspection of industrial parts[END_REF] focus on the definition of scan path planning that optimizes the quality of the point cloud collected. A point of view is considered as optimal if the view angle, defining the scanning direction, and the scanning distance lie within admissible intervals. These intervals ensure the digitizing noise to be less than the expected threshold. Their approach is specific to laser-plane sensors. Rémy [START_REF] Rémy | Contribution à l'automatisation du processus d'acquisition de formes complexes à l'aide d'un capteur laser plan en vue de leur contrôle géométrique[END_REF] decomposes his method in 3 main steps: calculation of the part visibility or local visibility, calculation of the global visibility and finally calculation of the actual visibility considering both the laser beam and the CCD matrix visibility. The part surface is approximated thanks to a STL meshing whereas the space is modeled as a meshed sphere. Each vector normal of the meshed sphere defines a point of view. For each facet of the part surface S i , the local visibility is obtained for a point of view p j if the angle between the directions of the point of view n j , and the local normal vector to the facet n i is less than 60°. This value ensures the quality of the scanning. Martins et al. proposed a voxel-based approach. The CAD model is reduced into a coarse voxel-map allowing the calculation of the local normal vectors, and defining the space volume to calculate collision free trajectory. In function of the point of view accessibility, voxel can be classified as Surface, Empty or Inside [START_REF] Martins | Automated 3D surface scanning based on CAD model[END_REF]. The volume is divided into slices the width of which is given by the laser stripe. For each slice, a collision free scan path is generated. Consecutive paths are linked together according to a zig-zag strategy to define the final trajectory. Considering a NURBS surface Raffaeli et al. propose to subdivide the surface into portions that are included in the scanner field of view, and for which the variations of the normal vectors do not exceed a threshold. Each portion is afterwards sampled [START_REF] Raffaeli | Context dependent automatic view planning: the inspection of mechanical components[END_REF]. The next step of the process is point clustering. The points that can be seen from the same point of view are gathered according two conditions: the maximum of the distance between 2 points belonging to the same cluster must be less than the scanner field of view, and the maximum angle between normal vectors must remain less than 90°.

Most of the literature methods are generally well adapted for a given sensor. The method proposed in the paper is generic and suitable for any type of sensor. The approach adopts the sensor standpoint starting from the minimum number of points of view that can be defined in relation with the sensor field of view and the part dimensions. This leads to the part modeling as a voxel-map wherein the voxel size is given by the sensor field of view. The analysis thus consists in finding portions of the surface visible by the sensor point of view. A quality criterion is also added to the notion of surface visibility.

VOXEL-BASED PATH PLANNING

The originality of the method is that it is generic for different types of digitizing sensors. For this purpose, each sensor is modeled as cone with regard to its field of view as displayed in Fig. 2. The cone axis defines the direction of the point of view, its basis corresponds to the dimensions of the field of view, and its height is given by the optimal digitizing distance. As scan path planning aims at minimizing the number of points of view leading to the complete part digitizing, the interest of the proposed approach lies in initializing the search process by considering a reduced number of points of view allowing the entire digitizing of the bounding part volume. This reduced number of points of view is defined in relation with the size of the sensor Field of View (FOV) and the part dimensions (Fig. 3).

Fig. 3: Representation of the initial set of viewpoints

The initial directions are first selected along the X, Y, Z-axes of the space. This initial set of viewpoints is equivalent to the part modeling as a voxel-map for which the voxel size (L opt ) is equal to the size of the sensor FOV (Fig. 3). Visibility analysis is performed from this modeling. Before giving details of the approach, a series of definitions is required.

Normal vector to the voxel

A voxel is 3D pixel, most generally a cuboid. The representation of a digitized cloud of points as a voxel-map allows the identification of empty or surface voxels. In this paper, the concept of voxel-space introduced in [3] is adopted. Therefore it is possible to attach attributes, such as the normal or the barycenter, to each voxel in function of the surface portion included inside.

As the CAD model is known, it is easy to retrieve the STL representation of the surface. A surface voxel includes a portion of the object surface that means a discrete number of facets. The vector normal is defined considering the mean value of the normal vectors to each facet included in the voxel (Fig. 4):

1 1 . . N kk k v N kk k Sn n Sn = = = ∑ ∑ r r r (2.1)
where, v n r is the vector normal to the considred voxe, k n r is the vector normal to the facet k included within the voxel, k S is the facet area, and N is the nimber of facets included in the voxeL. As the vector normal will be sensitive to the number of facets included into the voxel, the normal vector must be assessed to check if it is representative or not of the surface portion that is included into the voxel. For this purpose, the vector normal consistency is checked as proposed in the next section.

Consistency analysis

The normal vector previously calculated is afterwards used to describe the surface portion included in the voxel. To ensure that this vector is representative of the surface in terms of visibility, a consistency analysis is performed. Actually, the surface portion may present some abrupt local curvature variations, that means normal variations non coherent with the mean value, and thus with the visibility evaluation (Fig. 4).

Fig. 4: Consistency analysis

The consistency analysis is performed considering the consistency between each local normal vector normal to the mean normal vector normal. Let us consider a facet k. If k n r represents its local normal vector, the angle k θ between the mean normal and the local normal can be defined as follows:

. cos .

kv k kv nn arc nn θ  =    rr rr (2.2)
Consistency is assessed if:

kt h r e s h o l d θθ ≤ (2.
3)

The value of threshold θ is defined from the cone used to describe the sensor. Usually, the value of threshold θ is chosen equal to a quarter of the angle defining this cone.

Voxel qualification

A voxel V will be said seen from a point of view (,) pMd r if it is the first voxel intercepted by the line issued from M and directed along d r (Fig. 3). Nevertheless, although a voxel is seen by a point of view, the digitizing quality cannot be reached. Indeed, numerous studies enhanced the strong influence of the view angle on the quality of the collected points. In fact, if the angle between the digitizing direction and the vector normal to the voxel is greater than a threshold, the quality of the collected points may be affected in terms of trueness or in terms of digitizing noise [START_REF] Mehdi-Souzani | Contact less laser-plane sensor assessment: toward a quality measurement, IDMME-Virtual Concept[END_REF] [START_REF] Audfray | Assistance to automatic digitizing system selection for 3d part inspection ASME[END_REF]. In this direction a qualification of the seen voxels considering quality criteria is proposed. This qualification is possible as the analysis is performed only for voxels for which the vector normal is consistent with the surface portion included.

METHOD DESCRIPTION: VOXEL2SCAN

The algorithm aims at defining the reduced set of points of view allowing the complete part surface digitizing with a given digitizing quality. The algorithm is initialized considering a set of initial points of view directed along the 3 principal axes of the part coordinate system, with the objective of scanning large portions of the part surface. The part is modeled as a voxel-map for which the voxel width is given by the size of the Field of View (FOV). Non digitized details are collected thanks to a refinement of the voxel-map and by considering additional points of view. The algorithm consists of 4 main steps:

1. Definition of the voxel-map of the surface and determination of the set of initial points of view 2. Adaptive refinement of the voxel-map to ensure consistency of the normal vectors 3. Qualification of the voxels (well-seen, poorly-seen, not-seen) according to the initial set of points of view 4. Determination of additional points of view for not-seen voxels It is important to note that a process of collision checking is done during steps 3 and 4, in order to obtain a final set of collision free points of view. The main idea here is to propose an algorithm, generic, fast and simple, that leads to a set of collision free points of view allowing the part surface digitizing within a given quality. The 4 steps are detailed in next sections.

Definition of the initial voxel-map and the initial set of points of view

The part surface to be digitized is modeled thanks to a STL meshing. This type of representation is simple enough, and allows the calculation of local properties, such as the normal vector as detailed in 2.1. The determination of the initial set of points of view is defined as follows:

-Voxelization of the part: the size of each voxel is given by the width L opt of the FOV as proposed in figure 2; each voxel includes a portion of the meshing that serves to calculate the normal vector according to Eqn. (2.1).

-Determination of the initial set of directions of view: the directions are chosen along the axes (X, Y, Z) that compound the coordinate system linked to the voxel-map, in both the negative and the positive directions; these 6 initial directions are chosen to simplify the first visibility analysis.

-Determination of the set of points of view: the initial voxel-map is included in a parallelepiped the faces of which are perpendicular to the 6 initial directions (Fig. 3); this parallelepiped is enlarged by a distance equal to the optimal digitizing distance D opt , and the center of each voxel is projected onto the nearest face of the enlarged parallelepiped defining the origin of the point of view M; the direction of the point of view d r is given by the vector normal to the pane.

At this stage the set of the initial points of view (M j , i d r

) is defined.

Adaptive refinement

For each voxel previously defined, the consistency analysis is carried out. If the consistency is not ensured, the voxel size is divided by two, and the consistency is assessed for each one of the two new voxels. This stage is performed iteratively until the consistency is reached. Nevertheless, a shutoff parameter must be considered here to avoid too small voxel-size. Therefore, the process ends if the voxel size is equivalent to the average mesh size of the STL surface representation, even though the consistency is not reached. Each new voxel generates a new point of view as defined in the previous section (by projecting its center onto the principal nearest plane). If the vector normal is not consistent with the surface portion when the shutoff parameter is reached, the voxel is identified as "Difficult" and an orientation is proposed in function of the voxel neighborhood.

It is important to notice here that the refinement is local and well-adapted to surface portions presenting small curvature radius variations.

Qualification of the voxels

From the set of initial points of view, two different voxel lists are created in function of each voxel qualification: a list of the seen voxels L S and a list of the L NS not-seen voxels. At the algorithm initialization the first list is empty while the second includes all the voxels. For each direction of view i d r , the approach is the following:

1. Determination of the voxels belonging to L NS that are visible from i d r

; these voxels are classified in L S . 2. Qualification of the seen voxels from according to the approach detailed in section 2.3; the voxels poorly-seen are removed from L S and classified in L NS . 3. For each voxel V j belonging to L S , the position of the point of view M j is the offset of the center of the voxel by a distance equal to D opt according to the considered direction; to the voxel V j corresponds the point of view p j (M j , i d r ) 4. Collision checking for the point of view p j ; a collision arises if a sphere, the radius of which is given by the sensor dimensions, centered at M j intersects at least one voxel; in this case, the voxel V j is removed from L S and classified in L NS 5. Selection of another direction

1 i d + r
and reiteration of the process from step 1 At this stage, to each voxel V j belonging to L S corresponds a point of view p j .

Additional points of view

At the end of the previous step, the list L NS may be not empty, that means that some voxels remains not-seen or poorly-seen. Finding additional directions is thus necessary. Various methods are possible. For instance, optimized points of view could be defined by considering for each voxel the best view with regard to its vector normal. However, this could add numerous points of view that could make the trajectory planning difficult. Therefore, additional directions of view simply result from the intersection of the bisector planes of the coordinate system with third one. This solution works as for most of the sensors, half the admissible view angle α adm /2 is greater than 45°. These directions are added for all the voxels include in L NS until the list becomes empty. In the case of surfaces containing important variations of local normal vectors, the proposed algorithm may lead to a great number of voxels with small dimensions due to the refinement. The dimensions may be less than the size of the FOV. This involves that the distance between two consecutive points of view may be less than the width of the FOV. In this case, a portion of the surface may be digitized several times. To avoid this problem, the initial positions calculated from the center of the voxels of the initial voxelization are kept for these additional directions.

Once, all the points of view are determined, the sensor trajectory must be defined. As this trajectory strongly depends on not only the sensor but on the whole digitizing system used, the construction of the trajectory will not be discussed here.

APPLICATION

To assess the method, the previous algorithm is applied to a test surface. The surface consists of a plane, a truncated cone bonding by a spherical portion and a ruled surface. This test surface has been designed as multiple points of view are required for its complete scanning. The method is applied to two different sensor technologies: a laser-scanner (Kreon Zephyr KZ 25) and a light-projection system (Atos -GOM) with the parameters displayed in Table 1. The choice of the admissible angle of view is given thanks to a first sensor assessment process to evaluate its quality in terms of trueness and of noise [START_REF] Zuquete-Guarato | Towards a new concept of in-line crankshaft balancing by contact less measurement: process for selecting the best digitizing system[END_REF]. Results are shown in figure 6. The refinement remains local due to the weak surface complexity. However, it is interesting to observe that the proposed method well operates regardless of the sensor technology used and provides a set of points of views ensuring the desired quality of the collected points. For both systems, the admissible angle of view is 60°. This value ensures a noise less than 3 µm for the ATOS system and around 15 µm for the KZ25 sensor. In these conditions, the trueness can be evaluated at 3.5 µm for both digitizing systems [START_REF] Zuquete-Guarato | Towards a new concept of in-line crankshaft balancing by contact less measurement: process for selecting the best digitizing system[END_REF]. As expected the number of points of view is reduced for the large FOV sensor. The obtained point density is relatively homogeneous as overlap zones are small (Fig. 7). The comparison of the measured points with the CAD model shows that the quality obtained is also homogeneous. The main deviations are located near corner radii lie and are likely due to the machining operation. However the digitized point cloud presents digitizing holes corresponding to non-measured areas of the surface. Actually, as the model of the sensor is reduced to one point and a direction, it cannot account for the visibility of each device of the sensor (visibility laser and visibility camera). These digitizing holes illustrate one limit of the method linked with the sensor model simplicity. Fig. 7: Result of the part digitizing with the KZ25 sensor Nevertheless, more than 95% of the surface is measured with an expected trueness around 3.5µm and an expected noise around 15 µm. The measured surface is smaller than the nominal one: 0.024 m 2 for the CAD model for 0.025 m 2 for the measured one. After registration of the measured data on the CAD model, the mean deviation is evaluated to 0.0007mm. Deviations are mainly due to machining defects near the corner radius. Some digitizing errors can be observed close to the external edge of the surface.

CONCLUSION

An original method to determine optimal scanning path is proposed based on the size of the sensor field of view. Starting from a reduced set of points of view admissible with regard to the field of view and the part dimensions, the algorithm checks the validity of the points of view in terms of visibility of the surface and digitizing quality. The visibility analysis relies on an adaptive voxel-map of the surface which limits significantly computational time. As illustrated through two different digitizing systems, the method applies for any type of sensors leading to the set of points of view allowing the part digitizing. The scan trajectory is afterwards defined in function of the digitizing system used. The application enhances the efficiency of the algorithm to find, with a minimum calculation effort, the set of points of view leading to the digitizing of the part surface with a given quality. The simplicity of the sensor modeling shows its limits as it does not account for the actual sensor visibility. An improvement of the method is actually in progress to account for this limit.

Fig. 1 :

 1 Fig. 1: Definition of a point of view

Fig. 2 :

 2 Fig. 2: Sensor modeling in function of the sensor Field of View (FOV)

4 )Fig. 5 :

 45 Fig. 5: Voxel qualification

1 .

 1 Algorithm parameters for the two sensors

Fig. 6 :

 6 Fig. 6: Points of view for the ATOS system (left) and the KZ25 sensor (right)