Vincent Guillet
email: vincent.guillet@univ-orleans.fr

Guy Lamarque
email: guy.lamarque@univ-orleans.fr

Philippe Ravier
email: philippe.ravier@univ-orleans.fr

Christophe Léger
email: christophe.leger@univ-orleans.fr

Reduction of the impulsive noise emitted by a resistive load commutation

Keywords: Commutation, impulsive noise, Signal-to-Noise Ratio (SNR) improvement

In power line or any cable communications, impulsive noise is assumed to be the most difficult noise to filter. In particular, it is impossible to predict when a non periodic asynchronous impulsive noise will appear. So, under such noise conditions, the OFDM 1 symbol generally used in PLC 2 , is corrupted. Signal-to-noise ratio is generally improved by detecting and/or filtering such noise to overcome this problem. Unfortunately, this leads to heavy detection and computing time in comparison with the disturbance duration. In this paper, we propose a new approach that consists in controlling the instant of commutation of a load. This approach allows reducing up to 15 dB the asynchronous impulsive noise emitted by a radiator in a 75-μs typical OFDM symbol duration.

I. INTRODUCTION

Nowadays, the great majority of systems needs to communicate without extra wire. When the systems are portable, it is quite common to use wireless communication facilities. However, when the systems require the power supply, the communication over power line appears to be a good alternative solution. Actually, this communication channel does not need any additional wire, which is effortless and cheap to set a network in existing buildings. Moreover, it offers Triple Play facilities (Internet, VoIP, Multimedia) thanks to a high PHY 3 bit rate, in the range on few hundreds Mbps (Mega bit per second) raw data, and can be achieved through around one hundred meter cable. Thus, achieving such a high bit rate needs allocating the authorized bandwidth in an optimal way. An appropriate modulation like OFDM 1 which has a good spectral efficiency allows taking into account the overall bandwidth. Nevertheless, it is very important to have both a good knowledge of the topology of the communication channel [START_REF] Banwell | A new approach to the modeling of the transfer function of the power line channel[END_REF] [START_REF] Anastasiadou | An experiment setup for characterizing the residential power grid variable behavior[END_REF], and the associated noises.

Considering the power line communication channel, it is well known that the noise does not agree with the classical Additive White Gaussian Noise (AWGN) model. The model below (Fig. 1), presented in [START_REF] Zimmerman | An analysis of the broadband noise scenario in powerline network[END_REF], is generally preferred. In this model, the communication channel is mainly disturbed by a set of five noise sources. The background noises [START_REF] Liu | Broadband characterization of indoor powerline channel[END_REF] (narrow band and colored noises) are frequency localized. Conversely, the impulsive noises [START_REF] Chan | Amplitude, width, and interarrival distributions for noise impulses on intrabuilding powerline communication networks[END_REF], [START_REF] Phillips | Development of a statistical model for power-line communication channels[END_REF] are time concentrated. Among the impulsive noise class, some are periodic (periodic asynchronous/synchronous impulsive noises) and the others are not (asynchronous impulse noises).

Among all these types of noise, the asynchronous impulsive one is probably the most difficult noise to deal with. Numerous methods allow to detect and/or filter such noises [START_REF] Chan | Amplitude, width, and interarrival distributions for noise impulses on intrabuilding powerline communication networks[END_REF] [7] [START_REF] Umehara | Modeling of impulse noise for indoor broadband power line communications[END_REF]. However, these methods are computing time consuming and not adapted for real time implementations [START_REF] Wouters | Real time implementation on fpga of an ofdm based wireless lan modem extended with adaptive loading[END_REF] [START_REF] Dick | Fpga implementation of an ofdm phy[END_REF].

As far as the topology is concerned, the main element to take into account is the model of the cable. The cable used for fix installations is rigid. Generally, rigid cable is used for fix installation because it is cheaper than flexible one, which is used in complement or for extensions for its installation convenience.

In this paper, a new approach based on the reduction of the aperiodic asynchronous impulsive noise emitted by a resistive load commutation is proposed. The signal-to-noise ratio in such a case is improved by controlling the instant of commutation of the load.

The rest of this paper is organized as follow. Section II introduces the commutation system and the way to analyze the noise emitted. Section III presents the results of the commutation. Section IV demonstrates that the response of the impulsive noise is not only due to the commutation instant but also depends on the topology of the network.

II. MEASUREMENT PROTOCOL

The study of the commutation effects of a load is carried out using the schematic presented in figure 2) and described in the following paragraph.

A. The cable

Cables installed in buildings do not use a global sheath for the two or three copper conductors. However, a 3G cable (See Table I) was used in order to fix the inter-conductors distance.

B. The load and the commutation

The load consists in a 240V / 500 W radiator placed at the end of a cable. This type of load is purely resistive. Its commutation is done using a relay technology. This technology is equivalent to the bilame used in a thermostat of a radiator. A digital dephasing system synchronous to the main is used to evaluate the importance of the instant of the commutation (Fig. 3). The zero-crossing circuit detects the instant of a zero amplitude value of the main power supply. It is carried out using an optocoupler [START_REF] Dostert | Power-line comunications[END_REF].

At this precise instant, a digital counter implemented in a CPLD 4 (Xilinx XC9572 [12]) starts counting. A 32.768 kHz system clock increments the counter value asserting a precision of t ≈ 30.5 μs on the commutation instant. When the delay value is reached, the desired state of the commutation takes effect.

C. Analysis

The analysis of the channel is carried out using a standard capacitive coupling circuit. The circuit is designed to reject low frequencies (50 Hz), that carry the main power supply. It lets high frequencies passing through. The bandwidth is in the range of 100 kHz to 35 MHz. The analyzer is a dual input acquisition board GAGE CompuScope 12400 [START_REF]Gage compuscope 12400 specifications[END_REF], with 12 bits of resolution, and sampling at a fixed rate of 200 MHz. Due to such data flow, the on-board memory is used preventing the loss of samples during the acquisition. A hundred times attenuator is used to avoid saturation of the acquisition tool. The input impedance of this attenuator is 50 Ω. Signal processing is carried out using MATLAB environment in differed time. The measurement includes an overvoltage protection using a sufficient bandwidth to insure a correct analysis.

In order to reduce the random part of commutation signal and the bounce of contact, ten acquisition sequences are carried out for each parameter (cable length and dephasing value).

III. COMMUTATION ANALYSIS

The time duration of commutation pulses is from 1 μs to 50 μs. In few cases, like in the one presented in figure 4, some bounce can be observed. These bounces are due to the relay technology that has been used for the commutation system. Nevertheless, such bounces can also be observed in a bimetallic commutation employed in the thermostat of a radiator.

A time-frequency analysis is carried out using the Short-Time Fourier Transform (STFT) (Fig. 4). The square modulus of the STFT also called spectrogram is estimated using 1024points rectangular windows (corresponding to an about 5 μs window duration). The 1024-points FFT spectra are evaluated on sliding windows with 95% overlap between the slices.

The time frequency diagram (Fig. 4) clearly demonstrates that during the full pulse duration time, all frequencies, including authorized PLC bandwidth, are perturbed. Moreover, the disturbance duration (from 2 μs to 20 μs) is of the same order than the OFDM symbol duration of PLC standard protocols (from 10 μs to 100 μs). This consideration establishes that the disturbance time and the detection/processing time are equivalent. So, detecting pulses is not efficient. This is why we propose a new approach that controls the instant of commutation, in order to reduce the noise.

A. Measurement for a 30-meters cable

Figure 5 presents the noise PSD (Power Spectrum Density) obtained using a welch periodogram over 2048 samples with a rectangular window and without any overlap in a 75-μs typical OFDM symbol duration. The temporal delay is settled to 0, 12.5 and 25 percent of the power supply period. It respectively corresponds to commute the load when the power supply sinusoidal amplitude is equal to 0%, 70% and 100% of the maximum voltage value. A 30-meters cable is used. Figure 5 shows that the noise floor increases as a function of the delay of the commutation. In order to easily estimate the SNR improvement, an approximation model of the background noise is used [START_REF] Liu | Broadband characterization of indoor powerline channel[END_REF].

N (f ) = 10 a-b.f c (W/M Hz) (1)
The fitting curves using MATLAB function are presented in figure 6. The a, b and c parameters are given in table II.

On figure 6, it is obvious that the noise power is included between the noise level given by the curves at 0% and 25% instant commutation. Controlling the instant of commutation permits to fix the noise floor to its lower value. Doing so, the SNR can be improved by about 15 dB in a 75-μs typical OFDM symbol duration.

In order to ensure the independence to the length of the cable, the experimentation is also realized with a 70-meters cable.

B. Measurement for a 70-meters cable

The PSD obtained for a 70-meters cable is presented in figure 7. Figure 8 presents the fitting results.

The PSD and the fitting curves are estimated in the same manner as in section III-A. Table III gives the parameters obtained for the fitting. The fitting curves (Fig. 6 and Fig. 8) obtained for a 30meters and a 70-meters cables show an equivalent behaviour for a given commutation instant. Consequently, controlling the commutation instant may effectively reduce the SNR of about 15 dB, independently of the cable length. Regarding the figures 5 and 7, we notice ripples. The frequency ripple seems to depend on the cable length. Indeed, the frequency seems to increase with the cable length.

The ripple might be related to the line impedance mismatch.

In order to verify more precisely this hypothesis, a cable analysis is done in the next section.

IV. CABLE ANALYSIS

For this study, the same analysis circuit is employed, but the load is replaced by a generator and the commutation switch is left closed (Fig. 9). The generator is an arbitrary waveform generator (GAGE CompuGen 4300 [START_REF]Gage compugen 4300 specifications[END_REF]) using a white noise sampled at the frequency of 300 MHz. The output of the generator matches the full scale input of a line driver (OPA2674) avoiding saturation. Another coupling circuit is added on the line driver output.

The results for a 30-meters, 50-meters and 70-meters cable are presented in figure 10.

The line impedance mismatch will produce reflexions over the line. These reflected signals will be added to the measured signal after a propagation time given by the equation (2).

p = l v p = l C √ r , using v p = C √ r ( 2
)
where l is the length of the cable, v p is the propagation speed, C ≈ 3.10 8 m/s is the light speed in vacuum and r is the permittivity of the dielectric.

During a propagation time t p , it is possible to observe a k number of a 1 f sinusoidal signal period.

t p = k * 1 f (3)
A maximum value is obtained in f k for each integer value of k (4).

f k = k t p (4)
More over, the frequency distance between two local maxima is ΔF [START_REF] Chan | Amplitude, width, and interarrival distributions for noise impulses on intrabuilding powerline communication networks[END_REF].

ΔF = f k+1 -f k = k + 1 t p - k t p = 1 t p (5)
The measurement permits to evaluate the propagation time for a 50-meters cable and the dielectric cable permitivity using equation ( 2) giving equation [START_REF] Phillips | Development of a statistical model for power-line communication channels[END_REF] .

The frequency ripple value (the distance between two local maxima, ΔF ) for different cable lengths (0, 20, 30, 50, 70, 80 meters) are measured using zooms like the ones depicted on figure 11. The result is reported in figure 12, where equation ( 2) is also plotted. Both curves are superimposed demonstrating that the ripples are really due to the impedance mismatch. The measure of the frequency ripple value in figure 5 and figure 7 also verify the equation [START_REF] Anastasiadou | An experiment setup for characterizing the residential power grid variable behavior[END_REF].

Fig. 1 .

 1 Fig. 1. Power line channel noise description.

Fig. 2 .

 2 Fig. 2. Analysis of the commutation effects.

Fig. 3 .

 3 Fig. 3. Commutation tool.

Fig. 4 .

 4 Fig. 4. Time frequency analysis of a commutation.

Fig. 5 .

 5 Fig. 5. Comparison of the PSD waveforms as a function of three commutation instants for a 30-meters cable.

Fig. 6 .

 6 Fig.6. Fitting curves with respect to the background noise model for a 30-meters cable.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Comparison of the PSD waveforms as a function of three commutation instants for a 70-meters cable.

Fig. 9 .

 9 Fig. 9. Method to analyze the bandwidth of the cable.

Fig. 10 .

 10 Fig. 10. Ripples in the bandwidth.

  t

  Zoom for a 30-meters cable.

  Zoom for a 70-meters cable.

Fig. 11 .

 11 Fig. 11. Zoom on ripples in the bandwidth.

Fig. 12 .

 12 Fig. 12. Ripple ratio as a function of the cable length.

TABLE I CABLE

 I USED FOR THE TEST.

	Conductor type	Solid
	Manufacturer	La Triveneta Cavi
	Reference	U1000 R2V
	Class	1, Solid
	Insulation ( inner / outer )	PR / PVC
	Approximated permittivity 3.18 / 2.25 (Farad/meter) (inner / outer)

TABLE II FITTING

 II CURVES PARAMETERS FOR A 30-METERS CABLE.

	Commutation Switched instant Amplitude	a	b	c
	0%	0%	12.1 11.8 0.2
	12.5%	70%	22.3	9	0.27
	25%	100%	23.9	6	0.35

TABLE III FITTING

 III CURVES PARAMETERS FOR A 70-METERS CABLE.

	Commutation Switched instant Amplitude	a	b	c
	0%	0%	11.7 13.1 0.19
	12.5%	70%	22.4 9.4 0.27
	25%	100%	22.7 8.4 0.28

OFDM: Orthogonal Frequency Division Multiplexing.

PLC: Power Line Communication.

PHY: Lowest layer in the ISO model network communication.

CPLD: Complex Programmable Logic Device.

ACKNOWLEDGEMENT

This work has been developed in the frame of the French competitiveness cluster Science and Systems of Electrical (S 2 E 2 ) -SSAME GEDELOC -under the contract n2005-4020-06A-066.

CONCLUSION

In this study we have shown that the noise generated by asynchronous aperiodic impulsive sources can be easily and significantly reduced. Indeed, the SNR, in a 75-μs typical OFDM symbol duration, can be improved from up to 15 dB just by controlling the commutation instant of a resistive load.

Moreover, the asynchronous noise becomes synchronous. The pulses are always located at the zero crossing instant. Therefore, the bit loading can be optimized when the amplitude of the noise remains high. In that case, the data rate is increased, and the synchronization error is less corrupted. This work will be completed soon in both ways: trying different resistive and inductive loads and different network topologies.


Pierre Maurel
email: pmaurel@ece.umn.edu

Guillermo Sapiro

Facundo Mémoli

Alvaro Pardo

Renaud Keriven

Olivier Faugeras

Dynamic Shapes Average *

A framework for computing shape statistics in general, and average in particular, for dynamic shapes is introduced in this paper. Given a metric d(•, •) on the set of static shapes, the empirical mean of N static shapes,C 1 , . . . , C N , is defined by arg min C

. The purpose of this paper is to extend this shape average work to the case of N dynamic shapes and to give an efficient algorithm to compute it. The key concept is to combine the static shape statistics approach with a time-alignment step. To align the time scale while performing the shape average we use dynamic time warping, adapted to deal with dynamic shapes. The proposed technique is independent of the particular choice of the shape metric d(•, •). We present the underlying concepts, a number of examples, and conclude with a variational formulation to address the dynamic shape average problem. We also demonstrate how to use these results for comparing different types of dynamics. Although only average is addressed in this paper, other shape statistics can be similarly obtained following the framework here proposed.

Introduction

Understanding shape and its basic empirical statistics is important both in recognition and analysis, with applications ranging from medicine to security to consumer photography. The basic metrics and statistics of static shapes have been the subject of numerous fundamental studies in recent years, see for example [START_REF] Bookstein | Size and shape spaces for landmark data in two dimensions[END_REF][START_REF] Charpiat | Approximations of shape metrics and application to shape warping and empirical shape statistics[END_REF][START_REF] Cootes | Active shape models: Their training and application[END_REF][START_REF] Idryden | Statistical Shape Analysis[END_REF][START_REF] Miller | Groups actions, homeomorphisms and matching: A general framework[END_REF][START_REF] Small | The Statistical Theory of Shape[END_REF][START_REF] Toga | Brain Warping[END_REF][START_REF] Trouvé | Diffeomorphic matching problems in one dimension: Designing and minimizing matching functionals[END_REF] and references therein. In particular, given a metric on the set of static shapes (distance between two samples), the empirical mean shape of N static shapes, as well as other basic statistics, can be defined and computed. These are then used for diverse shape studies, from the recognition of particular objects to the detection of abnormalities in medical data. The purpose of this work is to extend this to dynamic shapes. This is fundamental for studies such as those involving gait, behavior, growth patterns, and all problems involving motion, deformations, and time-varying shapes.

Given N dynamic shapes Γ 1 (t), . . . , Γ N (t) (t stands for the time parameter, see Figure 1), we want to find M (t), the empirical mean of these shapes. This basic computation will be used throughout this paper as an example of how to perform statistics on dynamic shapes. One idea could be to simply perform static average among Γ 1 (t 1 ), . . . , Γ N (t 1 ), for each time instance t 1 , process that is clearly not efficient for every king of data. Indeed, the initial shape series might not be time-aligned (e.g., due to different growth rates in medical applications and different motion speeds in gait analysis). The dynamic shapes need to be properly aligned before any kind of shape statistics technique is applied. 1 This is exactly the role of the dynamic time warping (DTW), see Figure 1. This process is commonly used in speech recognition in order to time-align speech patterns to account for differences in speaking rates across speakers. It has also been used by a number of authors for gait analysis, but limited to the 1D path obtained by the tracking of particular joints. In this work we propose to combine DTW with results on static shape analysis to compute basic statistics on dynamic shapes. The framework here proposed is independent of the particular choice of static shape metric. This work deals with discrete time instances, while the extension to a continuous framework is discussed in the conclusions section.

Static Shape Averaging

In order to compute the mean of N static shapes, a distance on the set of the shapes is necessary (for examples of such metrics see [START_REF] Charpiat | Approximations of shape metrics and application to shape warping and empirical shape statistics[END_REF][START_REF] Kendall | Shape manifolds, procrustean metrics and complex projective spaces[END_REF][START_REF] Small | The Statistical Theory of Shape[END_REF][START_REF] Soatto | Deformotion: Deforming motion, shape average and the joint registration and approximation of structures in images[END_REF] and references therein). Once the metric is given, the empirical mean shape can be defined:

Definition 1 Let d(•,
•) be a distance on the set of shapes and C 1 , . . . , C N , N static shapes. The empirical mean M (C 1 , . . . , C N ) is given by

M (C 1 , . . . , C N ) arg min C N i=1 d(C, C i ) 2
The goal of this paper is to present a natural and easy to compute extension of this definition for dynamic shapes. Before doing this, let us briefly recall the second fundamental component of our approach, dynamic time warping.

Dynamic Time Warping

Dynamic time warping (DTW) is principally used in speech recognition to time-align speech patterns in order to account for differences in speaking rates across speakers. A distance between two speech patterns can then be computed by this technique in order to be able to compare them (see [START_REF] Parsons | Voice and Speech Processing[END_REF]). DTW can be adapted to deal with other types of signals as done in this paper for shapes.

When the two signals (A, B) to be matched are defined as sampled time functions, A = a 1 , . . . , a L ; B = b 1 , . . . , b M , the basic problem in DTW is to find two time warping functions f and g such that

T t=1 d(a f (t) , b g(t) ) 2
is minimized (here d(•, •) stands for the function measuring the discrepancy between two samples).

Computing these warping functions can be viewed as the process of finding a minimum-cost path through the lattice of points (a i , b j ) (i,j)∈{1,...,L}×{1,...,M } , starting from (1, 1) and ending at (L, M ) (see Figure 2), 2 where the cost of a path is defined by:

D(f, g) T t=1 d(a f (t) , b g(t) ) 2
and f and g are subject to the following constraints:

1. f and g must be monotonic:

f (k) ≥ f (k -1) and g(k) ≥ g(k -1)
2. f and g must match the endpoints of A and B:

f (1) = g(1) = 1, f (T ) = L and g(T ) = M
3. f and g must not skip any points:

f (k) -f (k -1) ≤ 1 and g(k) -g(k -1) ≤ 1 4. A limit in the maximum amount of warp is fixed by |f (k)-g(k)| ≤ Q, Q being the given "window width"
In the example in Figure 2, the time warping functions are:

f : 1 → 1, 2 → 2, 3 → 3, 4 → 4 5 → 5, 6 → 6, 7 → 6, 8 → 6 g : 1 → 1, 2 → 2, 3 → 2, 4 → 2 5 → 3, 6 → 4, 7 → 5, 8 → 6
At first glance, it would seem as if D(f, g) would have to be evaluated for a prohibitively large number of possible paths. Fortunately, dynamic programming brings this problem under control by noting that the best path from (1, 1) to any given point is independent of what happens beyond that point. Hence, if we call D(i k , j k ) the total cost of the best path from (1, 1) to (i k , j k ), this is the cost of the point (i k , j k ) itself plus the cost of the cheapest path to it:

D(i k , j k ) = d(i k , j k ) 2 + min legal(i k-1 ,j k-1 ) D(i k-1 , j k-1 )
By the subscript "legal (i k-1 , j k-1 )" we mean the minimum over all permissible predecessors of (i k , j k ). By constraints 1 and 3 above, there are only three legal predecessors: (i k -1, j k ), (i k , j k -1) and (i k -1, j k -1). Therefore we need to consider only three possibilities per lattice point (this if further constrained by point 4 above). Dynamic programming for solving the DTW problem (finding f and g) then proceeds in incremental stages (see [START_REF] Parsons | Voice and Speech Processing[END_REF] for the complete algorithm), achieving an optimal time complexity of O(P Q) (P is the number of initial frames and Q the "window width" from constraint 4). It means that we need to compute d(•, •), the distance between two static shapes, only O(P Q) times.

Dynamic Shapes Averaging

With the basic concepts on the mean of static shapes and dynamic time warping, we are now ready to describe the framework for dynamic shape average.

Basic Idea

We first define a dynamic shape as a sequence of static shapes (represented by any possible characterization): Definition 2 Let S be a set of static shapes (using any existing representation). A dynamic shape Γ is an ordered sequence of static shapes (C 1 , . . . , C T ) ∈ S T (T ∈ N is the length of the dynamic shape).

Although the above definition is given for discrete times, it can be extended to continuous space.

The idea now is to combine dynamic time warping and static shape averaging:

Definition 3 Given d(•,
•), a distance on the set of static shapes, and Γ 1 (t 1 ), . . . , Γ N (t N ), N dynamic shapes of respective length T i (i.e. t i = 1, . . . , T i ), their empirical mean is defined as

for 1 ≤ t ≤ T : Γ(t) M (Γ 1 (f 1 (t)), . . . , Γ N (f N (t)))
where f 1 , . . . , f N are N time-warping functions given by

(f 1 , . . . , f N ) = arg min f1,...,f N T t=0 µ(Γ 1 (f 1 (t)), . . . , Γ N (f N (t))) with µ(C 1 , . . . , C N ) = 1≤i<j≤N d(C i , C j ),
and M (C 1 , . . . , C N ) is the mean of static shapes (see Def.1).

In words, we start by finding (via DTW) optimal timecorrespondences between static shapes and after that we compute the average of these static shapes per time instance. The warping is such that the metric is minimized.

This definition suggests to consider the "distance" between two dynamic shapes as follows (there is no triangle inequality here): Definition 4 Given two dynamic shapes Γ 1 and Γ 2 , their "distance" is given by

δ(Γ 1 , Γ 2 ) = 1 T T t=0 d(Γ 1 (f 1 (t)), Γ 2 (f 2 (t))),
where d(•, •) is the selected metric for static shapes and f 1 and f 2 are the optimal time warping functions.

This definition will be used later to compare human motions.

Basic Improvements

With the simple use of DTW, the mean shape's length T will be greater than or equal to the maximum of the individual time lengths {T 1 , . . . , T N }. Therefore, the dynamic mean shape will always be longer than the initial shapes (in Figure 2, T 1 = 6, T 2 = 6, T = 8). In order to correct this, we add jumps in the final path.

When N = 2, define, for i ∈ {1, 2},

E i {( Γ(t), Γ(t + 1)) | f i (t) = f i (t + 1)}
where Γ(•) is the dynamic mean shape (from Definition 3). E 1 is the set of vertical segments and E 2 the set of horizontal segments in the graph representing the final path.

In Figure 2, E 1 = {( Γ(6), Γ( 7)), ( Γ(7), Γ( 8))} and E 2 = {( Γ(2), Γ(3)), ( Γ(3), Γ(4))}. These segments are responsible for the increase of the final length. Indeed, we have the simple relation (T is the length of the mean shape without jumps) :

T = T 1 + |E 1 | = T 2 + |E 2 |
We opt to replace every second pair in E 1 by its static average , then we do the same for the pairs in E 2 (see Figure 3). Each replacement decreases the length by one. Therefore T', the length of the mean shape with the jumps, becomes:

T = T - |E 1 | 2 - |E 2 | 2 = T 1 + |E 1 | -|E 2 | 2 T = T 1 + T 2 -T 1 2 = T 2 + T 1 2
The length of the final mean shape is then the average of the length of the two initial shapes.

In the general case (N dynamic shapes), it is also intuitive that we would like the length of the final mean shape to equal the average of the lengths of the N initial dynamic shapes. Therefore we now generalize the pairing process described above. Define, for any A ⊂ {1, . . . , N }:

E A {( Γ(t), Γ(t + 1)) | f i (t) = f i (t + 1) ⇐⇒ i ∈ A}
and, for any i ∈ {1, . . . , N }:

A i {A ⊂ {1, . . . , N } | i ∈ A}
Then the following relations, where T is the length of the mean shape without jumps, hold for any i ∈ {1, . . . , N } :

T = T i + A∈Ai |E A | and T - 1 N N i=1 T i = 1 N A⊂{1,...,N } |A| • |E A |
The right hand term of the previous equality can be eliminated by the following process: For every subset A of {1, . . . , N }, we choose a number |A|•|E A | N of pairs belonging to E A . 3 Then we replace each pair by their static average.The length of the final mean shape is then the average of the length of the N initial shapes. Figure 4 shows the mean shape for simple initial shapes and for N = 3.

Examples

For our experiments,we represented a static shape C by its distance function ψ(x) = min y∈C x -y and we used the following simple metric on the set of static shapes:

d(C 1 , C 2 ) = Ω (ψ 1 (ω) -ψ 2 (ω)) 2 dω
where ψ i (x) is the distance function to the shape C i . For this distance, C, the average shape, is the zero level set of

ψ(x) = 1 2 (ψ 1 (x) + ψ 2 (x))
As mentioned in the introduction, the framework here introduced is independent of the particular choice of the static metric d(•, •), and we have selected this simple one for demonstration purposes only.

The segmentation of the input pictures is done by simple thresholding ( [START_REF] Shapiro | Computer Vision[END_REF][START_REF] Sahoo | A survey of thresholding techniques[END_REF]), and the distance function ψ i for each shape C i is computed with the fast marching method (for details, see [START_REF] Helmsen | Two new methods for simulating photolithography development in 3d[END_REF][START_REF] Sethian | A fast marching level-set method for monotonically advancing fronts[END_REF][START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF]). Figure 5 shows an example (one frame) of an initial dynamic shape.

In figures 6 and 8, we present a number of frames from two initial video clips (dynamic shapes), followed by sampled frames from the average dynamic shape computed without using DTW, and finally sampled frames from the average dynamic shape computed with our technique. In Figure 10, we present some frames from three initial video clips (three walking men), followed by sampled frames from the mean dynamic shape computed without using DTW, and finally with our technique.

Using Definition 4, we can compare different dynamics, such as running vs. walking men. As observed in the table below, this function is five times greater between one running men and one walking men than between two running or two walking men. δ(•, •)/10 3 We choose these pairs uniformly spread in time.

Conclusion

A novel framework for performing shape statistics in dynamic shapes was described in this paper. The basic idea is to combine shape alignment with previously developed ideas from static shape studies. The shape alignment is based on dynamic time warping. The framework is independent of the metric between static shapes.

A number of directions are suggested by the line of research here initiated. First of all, other more advanced static shape metrics need to be used, including those that incorporate landmarks, found to be fundamental for medical applications [START_REF] Toga | Brain Warping[END_REF]. Once these advanced metrics are incorporated into our framework, we can proceed with more exhaustive experimentation, including 3D dynamic shapes. Of particular interest are the analysis and recognition of gait and the study of growth in medical applications.

In this paper we limited ourselves to the case of discrete time. In the continuous case, a variational formulation to address the dynamic shape average problem can be formulated as

arg min Γ,fi T 0 1≤i≤N [d(Γ(t), Γ i (f i (t)) 2 + H(f i (t)]dt
where f i are the time warping functions, and H represents some constraints on them (such as continuity, monotonicity, acceleration, etc). To this we can add time domain landmarks (e.g., by splitting the domain). These topics are the subject of current efforts in our group.

  Figures 7 and 9 show the corresponding DTW graphs.

Figure 1 .

 1 Figure 1. A simple example showing the importance of time alignment when performing shape statistics. The first two rows show two dynamic shapes Γ 1 (t) and Γ 2 (t). The following two rows show their mean: a) Computed without DTW-alignment, b) Computed with enhanced DTW-alignment. We clearly observe the need for the DTW step.

Figure 2 .

 2 Figure 2. Dynamic time warping example.

Figure 3 .

 3 Figure 3. Example of the introduction of jumps in DTW for N = 2.

Figure 4 .

 4 Figure 4. Example of mean shape, with jumps, for simple initial shapes and N = 3.

Figure 5 .Figure 6 .

 56 Figure 5. From left to right: Initial image, segmented shape, and distance function

Figure 7 .

 7 Figure 7. Graph corresponding to the DTW for the running sequences.

Figure 8 .

 8 Figure 8. Same as Figure 6 for two hands in motion.

Figure 9 .

 9 Figure 9. Graph corresponding to the DTW for the hands sequences.

Figure 10 .

 10 Figure 10. Example of three walking men. The three dynamic shapes are given first, followed by the mean without DTW (fourth row), and finally the mean with DTW, our proposed technique (last row). Once again, note the significant improvement when the time-warping is added to the shape statistics process.

  6 walk 1 walk 2 run 1 run 2

	walk 1	0	1.3	5.6	6.3
	walk 2	1.3	0	5.2	6.7
	run 1	5.6	5.2	0	1.1
	run 2	6.3	6.7	1.1	0

The topic of time alignment appears also in video (see for example[START_REF] Caspi | Parametric sequence-to-sequence alignment[END_REF] and references in there). The goals and techniques used there are completely different from the ones here presented. The use of our proposed framework for video alignment is the subject of future research.

Note that we use a i and b i both to denote the time positions and their corresponding values, the distinction clearly provided by the context.

work while visiting the University of Minnesota. We thank Renaud Keriven and Olivier Faugeras for suggesting and encouraging this visit. and Paul Thompson for interesting discussions on shape statistics. This work was supported by a grant from the Office of Naval Research ONR-N00014-97-1-0509, the Presidential Early Career Award for Scientists and Engineers (PECASE), and a National Science Foundation CAREER Award.


Guillaume Charpiat
email: charpiat@di.ens.frpmaurel@di.ens.fr

Olivier Faugeras
email: faugeras@sophia.inria.fr

Renaud Keriven
email: keriven@certis.enpc.fr

Pierre Maurel

Équipe Odyssée

DISTANCE-BASED SHAPE STATISTICS

This article deals with statistics on sets of shapes. The approach is based on the Hausdorff distance between shapes. The choice of the Hausdorff distance between shapes is itself not fundamental since the same framework could be applied with another distance. We first define a smooth approximation of the Hausdorff distance and build non-supervised warpings between shapes by a gradient descent of the approximation. Local minima can be avoided by changing the scalar product in the tangent space of the shape being warped.When non-supervised warping fails, we present a way to guide the evolution with a small number of landmarks. Thanks to the warping fields, we can define the mean of a set of shapes and express statistics on them. Finally, we come back to the initial distance between shapes and use it to represent a set of shapes by a graph, which with the technic of graph Laplacian leads to a way of projecting shapes onto a low dimensional space.

SHAPES AND SHAPE METRICS

Fully defining the notion of shape is beyond the scope of this article in which we use a limited, i.e purely geometric, definition. In our context we define a shape Γ to be the boundary of a regular and bounded subset of R n . Since we are driven by image applications we also assume that all our shapes are contained in a hold-all regular open bounded subset of R n which we denote by Ω. We suppose Γ to be a smooth codimensionone manifold of R n , and denote by S the set of shapes. We refer the reader to [1] for a more rigorous and complete analysis.

Since we want to be able to compare shapes, a way to quantify the similarity between them must be defined. One of the broadly used distance between shapes is the Hausdorff distance:

d H (Γ 1 , Γ 2 ) = max sup x∈Γ1 d Γ2 (x), sup x∈Γ2 d Γ1 (x)
where d Γ is the distance function to the shape Γ: d Γ (x) = inf y∈Γ d(x, y) . Some other distances are based on signed distance functions. The signed distance function to a shape Γ, denoted by dΓ , is equal to d Γ outside Γ and equal to -d Γ inside Γ. A possible shape distance is then the norm of the Sobolev space, W 1,2 (Ω), of square integrable functions with square integrable derivatives:

d W 1,2(Γ 1 , Γ 2 ) 2 = dΓ1 -dΓ2 2 L 2 (Ω,R) + ∇ dΓ1 -∇ dΓ2 2 L 2 (Ω,R n ) .

VARIATIONAL SHAPE WARPING

We assume that we are given a function E : S × S → R + , the energy. This energy can be thought of as a measure of dissimilarity between two shapes. Warping a shape Γ 1 into another one Γ 2 can be stated as the minimization of the energy E(., Γ 2 ) starting from Γ 1 , i.e., finding a family of shapes {Γ(t), t ≥ 0} with Γ(0) = Γ 1 and Γ(t) following some gradient descent toward Γ 2 .

Shape gradient

In order to define the gradient of the energy functional, the first step is to compute its Gâteaux derivatives in all directions, i.e., for all admissible velocity fields v : Γ → R n . Let us denote by G Γ E(Γ, Γ 2 ), v the Gâteaux derivatives of the energy function E(Γ, Γ 2 ) with respect to the shape Γ and in the direction v:

G Γ E(Γ, Γ 2 ), v = lim ε→0 E(Γ + ε v, Γ 2 ) -E(Γ, Γ 2 )
ε .

We would like to pick the gradient as the direction of steepest descent of the energy. But, to be able to assess the steepness of the energy, the deformation space needs to be equipped with an inner product. We model the space of admissible deformations as an inner product space F, , F . Under certain regularity conditions, there exists a vector w ∈ F such that:

∀ v ∈ F, G Γ E(Γ, Γ 2 ), v = w| v F .
We call it the shape gradient of E relative to the inner product , F and we note it w = D (F, , F ) Γ E(Γ, Γ 2 ). Usually F is taken as the set L 2 (Γ, R n ) of the square integrable velocity fields on Γ, and , F its associated inner product:

f |g L 2 = Γ f (x) • g(x) dΓ(x) .
In that case, we will only denote the gradient by D Γ E(Γ, Γ 2 ).

Equipped with some shape gradient, we can define the warping of a shape Γ 1 into another one Γ 2 as finding the family Γ(t) solution of the following Partial Differential Equation:

Γ(0) = Γ 1 ∂Γ ∂t = -D (F, , F ) Γ E(Γ, Γ 2 ) .
Natural candidates for the energy function E are the distances presented in the previous section. The Hausdorff distance is not Gâteaux differentiable. Yet, this problem can be solved using smooth approximation of this distance, denoted by dH (Γ 1 , Γ 2 ), which presents the advantage of being differentiable. This approximation is obtained using the following property (see [START_REF] Charpiat | Approximations of shape metrics and application to shape warping and empirical shape statistics[END_REF] for more details): for any f integrable on the compact set Γ, we have

lim β→+∞ 1 |Γ| Γ f (x) β Γ(x) 1/β = sup x∈Γ f (x) .
A few steps of the result of warping a silhouette of a fish onto another one, using the approximation of the Hausdorff distance as the energy function, are shown in Fig. 1.

Fig. 1.

Hausdorff warping of a fish onto another.

Generalized gradient and spatially coherent flows

Although mathematically well justified, the warpings induced by E = dH or E = d W 1,2 are not always completely satisfying: the obtained deformations do not seem to be the one a human observer would have chosen. To cope with this, a way to favor rigid (translations and rotations) and scaling motions is introduced in [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF] . This approach consists in changing the inner product used in the definition of the gradient. We decompose the deformation space L 2 into a sum of several mutually orthogonal linear subspaces: the subspace T of translations, the subspace R of rotations around the centroid, the subspace S of scaling motions centered on the centroid and the orthogonal complement of these subspaces G:

L 2 = T ⊕ R ⊕ S ⊕ G Applying different penalty factors to the different types of motions, we define a new inner product related to this decomposition and these penalty factors. The gradient of the Energy function relative to this new inner product can be easily deduced from the standard L 2 -gradient D Γ E (see [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF] for more details). The resulting warping promotes spatially coherent motions, keeping the warping of small details for the end of the evolution. In Fig. 2, we show a three-dimensional warping example from a teddy bear to Hayao Miyazaki's character Totoro. We use here the d W 1,2 distance defined in Section 1. The L 2 gradient descent is unable to give satisfying results. The modified gradient descent favoring rigid plus scaling motions leads to better results. Actually, only global coherent motions are promoted by this new gradient. The symptom of "unnatural" warping persists in the worst cases, such as complex shapes or shapes related by an articulated motion.

LANDMARKS

Landmarks are then necessary in worst cases. Provided by the user (anatomical landmarks), or automatically extracted (geometric landmarks), we assume that we are given p pairs of corresponding points on the initial and on the target shapes,

{(x 1i , x 2i ) ∈ Γ 1 × Γ 2 , 1 ≤ i ≤ p}.
We would like to use the information given by these correspondences to guide the evolution. We do this by adding a landmark term to the energy:

E tot (Γ(t), Γ 2 ) = E(Γ(t), Γ 2 ) +
Interestingly, the H 1 gradient can be obtained from the L 2 gradient by solving an intrinsic heat equation with a data attachment term (see [START_REF] Maurel | Reconciling landmarks and level sets: geometric shape warping and matching using generalized gradients and correspondence-augmented implicit representations[END_REF] for more details):

D H 1 Γ E tot is solution of ∆ Γ u = u -D L 2
Γ E tot where ∆ Γ denotes the intrinsic Laplacian operator on the surface, often called the Laplace-Beltrami operator. The solution of this equation coincides with that of the following variational problem:

arg min u Γ |u(x)-D L 2 Γ E tot (x)| 2 dΓ(x)+ Γ |∇ Γ u(x)| 2 dΓ(x)
and the H 1 gradient is finally a smoothed version of the L 2 gradient. Fig. 3 shows the warping between two hand shapes.

The energy E = d W 1,2 yields an unnatural warping. With three landmarks and spatially coherent flows, a satisfying warping is recovered.

MEAN, COVARIANCE AND MODES OF VARIATION

We have now developed the tools for defining several concepts relevant to a theory of stochastic shapes as well as providing the means for their effective computation. The first obvious concept is that of the mean of a set of shapes.

Definition 1 Given Γ 1 , • • • , Γ N , N shapes, we define their mean as any shape Γ that minimizes the function µ :

S → R + defined by µ(Γ, Γ 1 , • • • , Γ N ) = 1 N i=1,••• ,N E(Γ, Γ i )
An algorithm for computing approximations to a mean readily follows from the previous section: start from an initial shape and solve the PDE:

Γ t = -∇µ(Γ, Γ 1 , • • • , Γ N ) n.
We show some examples in Fig. 4.

We can now define the covariance of N shapes and their modes of variation (see [START_REF] Charpiat | Approximations of shape metrics and application to shape warping and empirical shape statistics[END_REF] for justifications).

Definition 2 Given N shapes Γ i , we note β i the normal velocity fields of the gradients of the functions Γ → dH (Γ, Γ i ) and β = 1 N N i=1 β i their mean. The covariance operator of these N shapes for their mean Γ is the linear continuous operator of L 2 ( Γ) defined by:

Λ(β) = i=1,N < β, β i -β > Γ (β i -β) . Definition 3 Let Λ be the N × N symmetric semi positive definite matrix Λ defined by Λij =< β i -β, β j -β > Γ. Let p ≤ N be its rank, σ 2 1 ≥ σ 2 2 ≥ • • • ≥ σ 2
p > 0 its positive eigenvalues and u 1 , • • • , u N the corresponding eigenvectors. Let u ij be the ith coordinate of the vector u j and v j be defined by

v j = 1 σj N i=1 u ij (β i -β). The velocities v k , k = 1, • • • , p
can be interpreted as modes of variation of the shapes and the σ 2 k 's as variances for these modes. Looking at how the shape varies for the kth mode is equivalent to solving the PDEs Γ t = ±v k (Γ) n with Γ(0, .) = Γ(.).

Examples of theses modes for the cases of the fingers and of the ten fishes is shown in Fig. 5.

GRAPH LAPLACIAN

We can interpret the previous section as a method to associate to any shape from a given set some coordinates in a low-dimensional linear space. Indeed, for any choice of dimension m, we can keep from a set of shapes only its m first eigenmodes of deformation (e 1 i m ) (the most significant ones, i.e., the ones with highest eigenvalues), which form an orthogonal family in the tangent space of the mean shape Γ embedded with its inner product. Then to a shape Γ we associate the warping field e Γ = ∇ Γd 2 ( Γ, Γ) defined on Γ and to this warping field we associate the coordinates e Γ |e i of its projection onto the family (e i ), i.e., we have associated to any shape Γ the coordinates of its projection in a low-dimensional linear space.

In fact, we do not need the computation of any warping field to be able to build a low-dimensional map of a set of shapes: even if these fields must carry much information, it is possible to build such a map when given only the distances between shapes, from a purely static geometrical point of view, thanks to the graph Laplacian technic [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF].

Let us denote by (Γ 1 i n ) a set of n shapes and consider a distance d on this set, for example the Hausdorff distance between shapes. We fix a positive integer K and search for the K nearest neighbors N i 1 l K of each shape Γ i for the chosen distance d. We then define a symmetric weight matrix W by

W i,j = δ i,j e -d(Γ i ,Γ j ) 2 2σ 2
where

δ i,j = 1 if i ∈ N j or j ∈ N i 0 otherwise
and we have chosen for σ the mean distance between neighbors:

σ = i,j d(Γ i , Γ j ) δ i,j i,j δ i,j .
Then, let us consider the symmetric negative semi-definite matrix L = W -D where D i,j = i W i,j δ i,j : it is a discrete approximation of the Laplacian operator. Thus, as explained in [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF], its eigenvectors F k of highest (negative) non-zero eigenvalues are the best functions from the shapes Γ i to R that could be used as coordinate system of the set of shapes. Consequently, we obtain a natural map in R m where each shape Γ i is represented by a dot with coordinates (F 1 k m (Γ i )). Let us try this approach on an artificial dataset. We build a set of rectangles with same center and width but different lengths and orientations, so there are two natural parameters we would expect the algorithm to find. Rectangles are randomly chosen such that the distribution of their corners is the uniform law in the authorized area (orientation between -π 6 and + π 6 , and length between 2 and 4 times the width). Results vary depending on the distribution density and the value of K: the higher the density, the better the results. Fig. 6 has been computed for 700 rectangles.

Let us now study the more complicated case of some different classes in a same connected component. We consider a set of 111 fish from the same database as before. The resulting map for the two first coordinates (see Fig. 7) shows some clusters of fish families.

CONCLUSION

We have presented several developments of our shape work, a generalization of our shape gradient, the extension of our framework that includes landmarks, and the application of a technique for dimensionality reduction that seems to offer interesting avenues for future work.

[1] M.C. Delfour and J.-P. Zolésio, "Shape analysis via distance functions: Local theory," in Boundaries, interfaces

Fig. 2 .

 2 Fig. 2. 3D shape warping with the L 2 gradient descent (top) and with a modified gradient descent favoring rigid plus scaling motions.

Fig. 3 .

 3 Fig. 3. evolution with E = d W 1,2 plus coherent flows plus three provided landmarks.

Fig. 4 .

 4 Fig. 4. Examples of means of several curves: a square and a circle (left), two ellipses (middle) and two hands (right).

Fig. 5 .

 5 Fig.5. The first mode of variation in the finger case and the first two ones for ten sample shapes of fish.

Fig. 6 .

 6 Fig. 6. Map from the graph Laplacian method for a set of rectangles whose length and orientation have been chosen randomly (K = 15).

Fig. 7 .

 7 Fig. 7. Two first coordinates for a set of 111 fish from different classes. The elements from each family are got together into clusters (K = 25).

p i=1 d(x i (t), x 2i ) 2 .where x i (t) is the ith landmark on Γ(t). In particular, we havex i (0) = x 1i .Formally, this energy yields Dirac peaks in the expression of the gradient of the energy. This is indeed not a good candidate for a gradient descent. The solution is again to change the inner product which appears in the definition of the gradient. Let H 1 (Γ, R n ) be the Sobolev space of square integrable velocity fields with square integrable derivatives. We consider the canonical inner product of H 1 (Γ, R n ):f |g H 1 = Γ f (x)•g(x)dΓ(x) + Γ ∇ Γ f (x)•∇ Γ g(x)dΓ(x)where ∇ Γ f and ∇ Γ g are respectively the intrinsic derivatives on Γ.


Pierre Maurel
email: pierre.maurel@ens.fr

Renaud Keriven
email: keriven@certis.enpc.fr

Olivier Faugeras
email: olivier.faugeras@sophia.inria.fr

Reconciling Landmarks and Level Sets

Shape warping is a key problem in statistical shape analysis. This paper proposes a framework for geometric shape warping based on both shape distances and landmarks. Our method is compatible with implicit representations and a matching between shape surfaces is provided at no additional cost. It is, to our knowledge, the first time that landmarks and shape distances are reconciled in a pure geometric level set framework. The feasibility of the method is demonstrated with two-and three-dimensional examples. Combining shape distance and landmarks, our approach reveals to need only a small number of landmarks to obtain improvements on both warping and matching.

Introduction

Understanding shapes and their basic empirical statistics is a fascinating problem that has attracted the attention of many scientists for many years [START_REF] Matheron | Random Sets and Integral Geometry[END_REF][START_REF] Kendall | A survey of the statistical theory of shape[END_REF]. Warping one shape into another is one of the keys leading to statistical shape analysis [START_REF] Dryden | Statistical Shape Analysis[END_REF]. Roughly speaking, the warping problem consists in transforming an initial shape into a target one: the result is the family of the intermediate shapes. Slightly different, is the matching problem, where a correspondence between two given shapes has to be established.

Introduced as a way to cope with interface evolution simulation, the level set method [START_REF] Dervieux | A finite element method for the simulation of Rayleigh-Taylor instability[END_REF][START_REF] Osher | Fronts propagating with curvaturedependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF] is based on an implicit representation of surfaces. The natural choice for the implicit representation is often the signed distance function to the closed surface. Consequently, the emergence of shape statistics in the implicit framework is not surprising. The pioneering piece of work considered the distance function as the only object of analysis: warping, matching, or statistical analysis were directly performed on the distance functions [START_REF] Paragios | Matching distance functions: A shape-to-area variational approach for globalto-local registration[END_REF][START_REF] Soatto | DEFORMOTION, deforming motion, shape average and the joint registration and segmentation of images[END_REF]. As we will see in Section 2 and 3 the Hausdorff distance [START_REF] Delfour | Shape analysis via distance functions: Local theory[END_REF][START_REF] Charpiat | Approximations of shape metrics and application to shape warping and empirical shape statistics[END_REF] can also be used as shape similarity measure.

Though, in the case of complex shapes, providing corresponding landmarks on both the initial and the target shapes reveals to be inevitable The natural way to guide an evolu-tion with landmarks is to try to minimize the distance between the landmarks on the evolving shape and the corresponding ones on the target shape. Again, this yields an irregular motion. In this paper, we present a novel usage of the generalized gradients introduced in [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF] that turns this motion into a regular and well posed one.

Remarkably, two recents advances in the level set method make our shape evolution compatible with it: first, a way to simulate a Partial Differential Equation (PDE) embedded on a surface [START_REF] Bertalmio | Variational problems and PDE's on implicit surfaces[END_REF], second, a way to deal with surface evolutions involving non normal velocities (and, but this is related, to track the surface points along time) [START_REF] Pons | How to deal with point correspondences and tangential velocities in the level set framework[END_REF].

We first review some shape distances and their usage for the shape warping problem. Then, after introducing the generalized gradient proposed by [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF], we present our landmark-guided warping. The next Section discusses the level set implementation of our method. The final Section shows two-and three-dimensional results and comparisons.

Shapes and Shape Metrics

In our context we define a shape Γ to be the boundary of a regular and bounded subset of R n , contained in the image Ω. We denote by S the set of shapes. We refer the reader to [START_REF] Delfour | Shape analysis via distance functions: Local theory[END_REF] for a more rigorous and complete analysis.

In order to compare shapes, a way to quantify the similarity between them must be defined. One of the broadly used distance between shapes is the Hausdorff distance:

d H (Γ 1 , Γ 2 ) = max sup x∈Γ1 d Γ2 (x), sup x∈Γ2 d Γ1 (x)
where d Γ is the distance function to the shape Γ: y) . The signed distance function to a shape Γ, denoted by dΓ , is equal to d Γ outside Γ and equal to -d Γ inside Γ. An other possible shape distance is then the norm of the Sobolev space, W 1,2 (Ω), of square integrable functions with square integrable derivatives:

d Γ (x) = inf y∈Γ d(x,
d W 1,2(Γ 1 , Γ 2 ) 2 = dΓ1 -dΓ2 2 L 2 (Ω,R) + ∇ dΓ1 -∇ dΓ2 2 L 2 (Ω,R n ) .

Variational Shape Warping

In this Section, we review the initial work of [START_REF] Charpiat | Approximations of shape metrics and application to shape warping and empirical shape statistics[END_REF] and its extension [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF].

We assume that we are given a function E : S × S → R + , the energy. This energy can be thought of as a measure of dissimilarity between two shapes. Warping a shape Γ 1 into another one Γ 2 can be stated as the minimization of the energy E(., Γ 2 ) starting from Γ 1 , i.e. finding a family of shapes {Γ(t), t ≥ 0} with Γ(0) = Γ 1 and Γ(t) following some gradient descent toward Γ 2 .

Shape gradient

In order to define the gradient of the energy functional, the first step is to compute its Gâteaux derivatives in all directions, i.e. for all admissible velocity fields v : Γ → R n . Let us denote by G Γ E(Γ, Γ 2 ), v the Gâteaux derivatives of the energy function E(Γ, Γ 2 ) with respect to the shape Γ and in the direction v:

G Γ E(Γ, Γ 2 ), v = lim ε→0 E(Γ + ε v, Γ 2 ) -E(Γ, Γ 2 )
ε For an inner product , F equipping the deformation space F , there exists a vector w such that:

∀ v ∈ F, G Γ E(Γ, Γ 2 ), v = w| v F
We call it the shape gradient of E relative to the inner product , F and we note it

D (F, , F ) Γ E(Γ, Γ 2 ). Usually F is taken as the set L 2 (Γ, R n ) of
the square integrable velocity fields on Γ, and , F its associated inner product:

f |g L 2 = Γ f (x) • g(x) dΓ(x) .
In that case, we will only denote the gradient by

D Γ E(Γ, Γ 2 ).
Equipped with some shape gradient, we can define the warping of a shape Γ 1 into another one Γ 2 as finding the family Γ(t) solution of the following Partial Differential Equation:

Γ(0) = Γ 1 ∂Γ ∂t = -D (F, , F ) Γ E(Γ, Γ 2 ) (1)
Natural candidates for the energy function E are the distances presented in the previous Section. The Hausdorff distance is not Gâteaux differentiable. Yet, this problem can be solved using an smooth approximation of this distance, denoted by dH (Γ 1 , Γ 2 ), which presents the advantage of being differentiable (see [START_REF] Charpiat | Approximations of shape metrics and application to shape warping and empirical shape statistics[END_REF] for more details).

Generalized gradient and spatially coherent flows

Although mathematically well justified, the warpings induced by E = dH or E = d W 1,2 do not reveal to be completely satisfying: the obtained deformations do not seem to be the one a human observer would have chosen. To cope with this, the same authors introduced in [3] a way to favor rigid (translations and rotations) and scaling motions. Their approach consists in changing the inner product used in the definition of the gradient (see [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF] for more details). Actually, only global coherent motions are promoted by this new gradient. We will see that the symptom of "unnatural" warping persists in case of complex shapes or shapes related by an articulated motion (see the hands example on Fig. 2).

Landmarks-guided warping

The energy

Landmarks are then necessary in many cases. Provided by the user (anatomical landmarks), or automatically extracted (geometric landmarks), we assume that we are given p pairs of corresponding points on the initial and on the target shapes, {(x 1i , x 2i ) ∈ Γ 1 × Γ 2 , 1 ≤ i ≤ p}. We would like to use the information given by theses correspondences to guide the evolution given by equation [START_REF] Bertalmio | Variational problems and PDE's on implicit surfaces[END_REF]. We do this by adding a landmark term E L to the energy:

E tot (Γ, Γ 2 ) = E(Γ, Γ 2 ) + E L (Γ, Γ 2 )
During an evolution, while forward correspondences may not exist (formation of shocks), backward correspondences are guaranteed: each point of the evolving interface comes from one point at time 0 (see [START_REF] Pons | How to deal with point correspondences and tangential velocities in the level set framework[END_REF]).

We note ψ t : Γ(t) → Γ 1 the family of functions giving for each point x of Γ(t) the point ψ t (x) on Γ 1 from which x comes. Let γ i (t) = ψ -1 t ({x 1i }) be the subset of Γ(t) coming from x 1i . Equipped with this correspondence, we are now able to define a landmark-based energy as the sum, for each landmark of Γ 2 , of the squared distance between this point and the corresponding set γ i (t):

E L = i d(x 2i , γ i (t)) 2 (2)
with the convention that the distance to an empty set is zero. Note that some landmarks might disappear (shock) or become a continuous infinity of points (rarefaction). Actually, we conjecture that, for reasonable choices of the landmarks, rarefaction does not happen with smooth curves. Yet, depending on the initial energy E, there might be some shocks, even with smooth curves.

In the sequel, we will suppose that either an initial landmark x 1i remains one point x i (t) (γ i (t) = {x i (t)}), or it disappears (γ i (t) = ∅). Under these hypothesis, the energy can be rewritten in the more classical way:

E L = {i,γi(t) =∅} d(x i (t), x 2i ) 2 (3)
keeping in mind that point x i (t) come from the backward correspondences ψ t .

Adapted gradient

Formally, the energy given by equation ( 3) yields Dirac peaks in the expression of the gradient of the energy:

D L 2 Γ E tot (x) = D L 2 Γ E(x)+ {i,γi(t) =∅} δ xi(t) (x)(x i (t)-x 2i )
where δ x denotes the Dirac function centered at point x. This is indeed not a good candidate for a gradient descent.

The solution here is inspired by [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF]. We change the inner product which appears in the definition of the gradient. Let H1 (Γ, R n ) be the Sobolev space of square integrable velocity fields with square integrable derivatives. We consider the canonical inner product of H 1 (Γ, R n ):

f |g H 1 = Γ f (x)•g(x)dΓ(x)+ Γ ∇ Γ f (x)•∇ Γ g(x)dΓ(x)
where ∇ Γ f and ∇ Γ g are respectively the intrinsic derivatives on Γ.

Interestingly, the H 1 gradient can be obtained from the L 2 gradient by solving an intrinsic heat equation with a data attachment term (see [START_REF] Maurel | Reconciling landmarks and level sets[END_REF] for more details):

D H 1 Γ E tot is solution of ∆ Γ u = u -D L 2 Γ E tot ( 4
)
where ∆ Γ denotes the intrinsic Laplacian operator on the surface, often called the Laplace-Beltrami operator. The solution of this equation coincides with:

arg min u Γ |u(x)-D L 2 Γ E tot (x)| 2 dΓ(x)+ Γ |∇ Γ u(x)| 2 dΓ(x)
(5) and the H 1 gradient is finally a smoothed version of the L 2 gradient.

Matching

Let us suppose that the warping process of Γ 1 into Γ 2 has converged. More precisely, we suppose there exists some time T such that Γ(T ) is very close to Γ 2 (e.g. E tot (Γ(T ), Γ 2 ) < ), and a way to assimilate points of Γ 2 to points of Γ(T ) (e.g. taking the closest point 1 ) . Then, the backward correspondence ψ T supplies a natural matching from Γ 2 to Γ 1 . This matching is not one to one if some points of Γ 1 have disappeared during the evolution (shocks).

Level set implementation

There is no need to introduce the broadly known level set method [START_REF] Dervieux | A finite element method for the simulation of Rayleigh-Taylor instability[END_REF][START_REF] Osher | Fronts propagating with curvaturedependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]. However, implementing our scheme in that framework requires two adaptations of the original method: implementing a PDE on an implicit surface and being able to track points during the evolution.

H 1 gradient

The H 1 gradient, solution of (4), is obtained from an iterative minimization induced by [START_REF] Dervieux | A finite element method for the simulation of Rayleigh-Taylor instability[END_REF]. Since the work introduced in [START_REF] Bertalmio | Variational problems and PDE's on implicit surfaces[END_REF], implementing a PDE on a surface is affordable in the implicit framework. The only hard point in our case could be the Dirac peaks in the data term. We indeed use a smooth approximation of them.

It should also be mentioned that, in the two dimensional case, the explicit solution of the equation ∆ Γ u = uv is known (see [START_REF] Maurel | Reconciling landmarks and level sets[END_REF] for more details) and can be used to avoid the iterative minimization giving u.

Point Correspondences

Because it codes interfaces with implicit representations, the original level set method can not follow the evolution of each point of the initial interface. Only the geometric location of the whole interface is recovered. Then, considered velocities are usually normal to the interface.

In our case, we need to follow the landmark points through the backward correspondences ψ t and to cope with the non normal velocity -D H 1 Γ E tot . We used a method proposed in [START_REF] Pons | How to deal with point correspondences and tangential velocities in the level set framework[END_REF] which enable to maintain an explicit backward correspondence from the evolving surface to the initial one.

Experiments

In the experiments showed here, we used the original energy E = d W 1,2 and test how our landmark-guided force modifies the warping and the final matching. Fig. 1 shows the warping of a rectangle into another one. The different parts of the curves are shown with different colors, so that their respective evolution can be followed. The initial warping without any landmark seems natural but it fails discovering the matching between the edges of the rectangles, a matching indeed recovered when providing landmarks. Fig. 2 shows the warping between two hand shapes. The energy E = d W 1,2 yields an unnatural warping. Adding spatially coherent flows makes the warping a bit better but still fails in some parts. With three landmarks only, both a satisfying warping and a good matching are recovered. Fig. 3 shows the warping of a teddy bear into a cartoon character. Without any landmarks, the top row evolution fails matching the ears and arms of the characters. The bottom row shows the evolution with four landmarks. Red spots allow to check a good matching between landmarks.

Conclusion

We propose a framework for shape warping based on both shape distances and landmarks. Our method is purely geometric and no extrinsic quantity like a space diffeomorphism has to be considered. Thanks to recent advances in the level set techniques, a level set implementation is possible, reconciling landmarks and the level set methods. Moreover, a matching between shapes is provided at no additional cost. Two-and three-dimensional examples, combining shape distance and landmarks, demonstrate the improvement brought by our approach on both warping and matching, even with a small number of landmarks. Further work includes investigating for a one-to-one matching between shapes, and a way to cope with other landmarks, such as curves on surfaces in R 3 .

Figure 1 .

 1 Figure 1. Warping of a rectangle shape into another one. Top row: evolution with E = d W 1,2 . Bottom row: same energy, with four provided landmarks, marked by color spots.

Figure 2 .

 2 Figure 2. Warping of a hand shape into another one. Top row: evolution with E = d W 1,2 . Middle row: same energy + spatially coherent flows. Bottom row: same energy + spatially coherent flows + three provided landmarks.

Figure 3 .

 3 Figure 3. Warping of a teddy bear into a cartoon character. Top row: evolution with E = d W 1,2 . Bottom row, first image: four landmarks (in blue) provided on the two shapes. Bottom row, remaining images: evolution with E = d W 1,2 plus the provided landmarks.

This could be a problem if the evolution gets stuck into some local minimum. Yet, we have never experienced this case.