Y Yang

S Chevallier

J Wiart

I Bloch

Automatic selection of the number of spatial filters for motor-imagery BCI

Common spatial pattern (CSP) is widely used for constructing spatial filters to extract features for motor-imagery-based BCI. One main parameter in CSP-based classification is the number of spatial filters used. An automatic method relying on Rayleigh quotient is presented to estimate its optimal value for each subject. Based on an existing dataset, we validate the contribution of the proposed method through a study of the effect of this parameter on the classification performance. The evaluation on testing data shows that the estimated subject-specific optimal values yield better performances than the recommended value in the literature.

Introduction

Neuro-electrophysiologic studies reveal that both real movement and motor imagery of a specific body part induce an electro-encephalogram (EEG) rhythmic attenuation termed event-related desynchronization (ERD) in the µ (8-13Hz) and β (13-30Hz) bands over corresponding functional regions in the sensorimotor cortex [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF]. Thus, the essential task of a motor-imagery based brain-computer interface (BCI) is to distinguish different spatial localizations of ERD for predicting different motor intentions. The common spatial pattern (CSP) algorithm is very effective in constructing optimal spatial filters that extract discriminative activity (i.e. ERD) and reduce feature dimensions in motor-imagery BCI [START_REF] Müller-Gerking | Designing optimal spatial filters for single-trial EEG classification in a movement task[END_REF]. This algorithm was firstly proposed for a binary discrimination and then extended to multi-class problems through various approaches (for details, see [START_REF] Dornhege | Increase information transfer rates in BCI by CSP extension to multi-class[END_REF]).

One main parameter of CSP-based classification is the number of paired spatial filters, which determines the features used in classification and therefore affects the classification result. Most researchers choose the value of this parameter just based on their experience and often use a constant value for all subjects, which ignores the potential individual differences. Although it was mentioned in [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis[END_REF] that this parameter can be alternatively determined via cross validation, this work neither provided any detail nor experimental validation. Moreover, using exhaustive searching strategy to find the optimal value of this parameter in the whole range of values increases the computational time, particularly when the dimension of data is very large. Thus, the method proposed in this paper includes two steps: 1) a criterion based on Rayleigh quotient is applied for pre-selecting the range of this parameter, 2) an algorithm based on cross validation is then employed for more precise estimation of the optimal value of this parameter in the pre-selected range. Based on an existing dataset, we validate the importance of the estimation through studying the effect of this parameter on the classification, and then verify the effectiveness of the proposed method by comparing the classification results using the estimated optimal values with those obtained using the recommended fixed value in existing work in both binary-class and multi-class problems.

Methods

Pre-selection of paired spatial filters

CSP is a data-driven approach to construct spatial filters, W = [w 1 , ..., w N ], which decomposes the N -channel EEG X = [x 1 , x 2 , ..., x N ] T into N uncorrelated filtered signals Z = [z 1 , z 2 , ..., z N ] T through the transformation z j = w T j X, (j = 1, 2, ..., N ) where w j is a generalized eigenvector that satisfies:

w T j C L w j = λ L j , w T j C R w j = λ R j , λ L j + λ R j = 1 (1)
where C L , C R ∈ R N ×N are the estimated covariance matrices of two classes (i.e. 'L' and 'R') of N -channel EEG signals, respectively. As λ L j + λ R j = 1, w j tends to yield a large variance of signal for one class and a small variance of signal for the other class. These contrary effects of w j on two classes contribute to the discrimination. Usually, w i and w N -i+1 according to i-th largest λ L and λ R , respectively, are used together as the i-th paired filters in CSP-based classification [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis[END_REF].

The discriminative activity S d and common activity S c between two classes are defined as

S d = C L -C R , S c = C L + C R ,
respectively. Thus, the ratio between discriminative activity and common activity projected on the w j spatial filter is the (Rayleigh quotient ) R(w j ) [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis[END_REF] and is obtained by:

R(w j ) = w T j S d w j /w T j S c w j = λ L j -λ R j /(λ L j + λ R j ) = 2λ L j -1 (2)
For the i-th paired filters, F D(i) = R(w i ) + R(w N -i ) reflects their effectiveness in extracting the discriminative components from the original signal [START_REF] Lv | Common spatial pattern and particle swarm optimization for channel selection in BCI[END_REF]. Usually the first m pairs of spatial filters according to the m largest F D(i) are used. Too small or too large values of m will lead to poor classification performances (see Sect. 3.2), so that the optimal value of m should be estimated for each subject.

A too small F D(i) (typically F D(i) < 0.1) indicates that the i-th paired filters have a very weak ability of extracting discriminative components, and cannot improve classification results (see Sect. 3.2). As all paired filters are sorted in descending order of F D(i), the F D(i) values are used as a pre-selection criterion to shrink the range for seeking the optimal m value.

Refined estimation of the optimal number of paired filters

The optimality criterion for selecting paired spatial filters must satisfy two properties: (1) the number of paired spatial filters must be minimal, (2) it must yield the classification result that is equal or comparable to the best one, i.e. such that there is no statistical difference between them or their difference is less than a tolerance δ (δ = 0.015 in this paper). Here, the classification performances are evaluated via the kappa coefficient : κ = (P o -P e )/(1 -P e ), where P o is the observed agreement between classifier and dataset labels, and P e is the chance level for agreement (i.e. P e = 0.5 for binary-class problems, P e = 0.25 for four-class problems). Thus, a larger κ value indicates a good classification result [START_REF] Dornhege | Toward brain-computer interfacing[END_REF].

Assuming the number of paired spatial filters with F D(i) ≥ 0.1 is M , the optimal number of paired spatial filters is evaluated by checking each possible m value (m ≤ M ) to see whether its corresponding κ value is significantly larger than others obtained for smaller values of m. The paired difference test (paired t -test) is employed for the significance analysis [START_REF] Woolson | Statistical methods for the analysis of biomedical data[END_REF]. In this case, if several m values yield equal or comparable classification results, the smallest one will be chosen as the most optimal. The algorithm of this procedure is described below. The optimal parameter m opt is estimated offline from the training data for each subject, and then applied to the testing data or on-line applications for the same subject. This strategy can be extended to multi-class problems (see Sect. 3.4).

Experimental validation

Data description

The data used in this work are from BCI competition IV dataset IIa [START_REF] Brunner | Spatial filtering and selection of optimized components in four class motor imagery EEG data using independent components analysis[END_REF], which contains one training session and one testing session of 22-channel EEG data from 9 subjects who performed four classes cue-driven motor imagery (left hand, right hand, both feet and tongue). Details about this dataset can be found in the associated technical document 1 . In this paper, we first use the data of left and right hands to investigate the effect of m on classification in Sect. 3.2 and then test the proposed automatic estimation strategy on binary-class (left vs. right hands) in Sect. 3.3 and multi-class data (the full dataset) in Sect. 3.4.

Sensitivity analysis of the number of paired spatial filters

A broad frequency band of 8-30Hz (µ and β bands) and the segment of 0.5-2.5s of EEG data after the cue onset were used in this study for calculating the transformation matrix W in CSP, and F D(i) value for each paired spatial filters, and for training the classifier [START_REF] Müller-Gerking | Designing optimal spatial filters for single-trial EEG classification in a movement task[END_REF]. The Fisher's linear discriminant analysis (LDA), which is classically used with CSP, was employed here for the classification [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis[END_REF]. The effect of the number of spatial filters was studied on the training data using 100 repetitions of 10-fold cross-validations. The classification performances were measured by κ value. Algorithms of CSP, classifier training and evaluation (including calculating κ value) are performed with the BioSig toolbox2 .

The effects of the parameter m on the classification results and F D(i) value of each paired spatial filters for all subjects are shown in Figure 1. It can be observed that (1) the performance of CSP-based classification is not proportional to m but has significant variations depending on m for all subjects; (2) the sensitivities of classification results to m are different between subjects: some (i.e. subjects 8, 9) are relatively low but most are relatively high; (3) adding the paired spatial filters with F D(i) < 0.1 does not improve the classification results: e.g. for subject 1, M = 6 and the κ value decreases if m > 6 is used; for subject 8, M = 7 and the κ value remains stable when m > 7. Those results prove that (1) it is critical to choose a right m value for each individual in CSPbased classification; (2) it is reasonable to estimate the optimal m in the range of [1, M ], where M is the number of the paired spatial filters with F D(i) ≥ 0.1.

Comparison for binary-class discrimination

Table 1 lists the estimated m opt values learned from the 100 repetitions of 10fold cross-validations in the training data of two classes (i.e. left hand, right hand) and provides a comparison of the evaluation results on the independent testing data using m opt and the classical value (m = 3) recommended in [START_REF] Müller-Gerking | Designing optimal spatial filters for single-trial EEG classification in a movement task[END_REF][START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis[END_REF]. The m opt value varies for different subjects and the classification performances are better than those with the recommended value. For subjects 8 and 9, whose sensitivities to m are relatively low, one pair of filters can already yield fine performance, while others may need more pairs of filters. Table 1: Estimated m opt values were obtained by 100 repetitions 10-fold crossvalidations on the training data. The computational time was always less than 32s on a 2.66GHz PC with Matlab (2010Ra). The evaluation results were computed on the independent testing data using m opt and the recommended value (m = 3).

Extension to multi-class problem

One Versus the Rest (OVR) CSP is a multi-class CSP approach that computes W for each class against all others and then projects the EEG signals on all the 2m × P chosen spatial filters (P is the number of classes, here P = 4) to extract the features, and then performs a multi-class LDA classification [START_REF] Dornhege | Increase information transfer rates in BCI by CSP extension to multi-class[END_REF]. Based on the pre-selection procedure in Sect. 2.1, each W generates a M value, thus P × M values are obtained. The largest M value (M max ) is chosen as the upper limit of possible m opt . Then m opt is estimated based on the classification results in the range of [1, M max ] using Algorithm A and then applied to the independent testing data. The comparison of results obtained with m opt and with fixed recommended m is shown in Table 2 for the four-class problem of BCI competition IV dataset IIa. Using m opt leads to better performances than using the fixed recommended m. As we used the broad frequency band (8-30Hz) of EEG signal in this work, it is difficult to make a comparison with the 1 st placed winner in BCI competition IV who extracted features from multiple narrow bands and reported the results based on searching the largest Kappa over the entire time range of the testing data using a 2-s sliding window [START_REF] Chin | Multi-class filter bank common spatial pattern for four-class motor imagery BCI[END_REF]. However, it makes more sense to compare with the 2 nd placed winner 3 who used the same frequency band, in order to validate the interest of using subject-specific m opt with OVR CSP. The comparison showed that m opt with OVR approach in CSP-classification needs less classifiers (only one multi-class LDA) and generates better mean performance. Table 2: Estimated m opt (obtained in less than 90s) and independent evaluation in a four-class problem using m opt and fixed m, and comparison with the 2 nd placed winner in BCI competition IV who also used 8-30Hz data but fixed m (m = 4) and pair-wise approach with three LDA and one Bayesian classifiers.

Conclusion

The number of spatial filters used in feature extraction affects the classification results. An automatic strategy based on Rayleigh quotient and cross validation is proposed to estimate the subject-specific optimal m value. Experimental results show that the estimated optimal m values vary for different subjects and often yield better results than those obtained with the fixed recommended value for both binary-class and multi-class problems. The proposed strategy can be applied on the training data to estimate the optimal value of m for each subject and then use it for the long term on-line classification of the given classes for the same subject to achieve the best results.

Algorithm A :

 : Selection of the optimal number of paired spatial filters Let M denote the number of paired spatial filters with F D(i) ≥ 0.1 and let m ≤ M ; κ(m) is a set of κ for a given m evaluated with a 100 repetitions of 10-fold cross-validation (κ(m) ∈ R 1×100 ), κ(m) is the mean value (over the 100 components), t(a, b) represents the p-values of paired t-test between vectors a and b 1:mi ← 1; mj ← 2 2:while mj ≤ M do 3: if κ(mj ) > κ(mi) + δ and t(κ(mi), κ(mj)) < 0.05 then 4: mi ← mj 5: endif 6: mj ← mj + 1 7:endwhile 8:mopt ← mi 9:return the optimal parameter, mopt

Fig. 1 :

 1 Fig. 1: Effect of parameter m on the left vs right hand classification and F D(i) of each paired spatial filters for all subjects in the BCI competition IV dataset IIa. The horizontal line on the right plot indicates F D = 0.1.

  ) 0.75 0.22 0.96 0.40 0.11 0.35 0.70 0.94 0.86 0.59 κ (m = 3) 0.67 0.13 0.96 0.40 0.09 0.25 0.69 0.93 0.82 0.55

3

  http://www.bbci.de/competition/iv/results/index.html ) 0.72 0.30 0.74 0.48 0.20 0.25 0.75 0.71 0.68 0.53 κ (m = 3) 0.69 0.30 0.71 0.47 0.20 0.25 0.74 0.71 0.50 0.51 κ (2 nd ) 0.69 0.34 0.71 0.44 0.16 0.21 0.66 0.73 0.69 0.52

http://bbci.de/competition/iv/desc_2a.pdf

http://biosig.sourceforge.net/

* This work was partially supported by grants from China Scholarship Council and Orange Labs.


Yuan Yang

Sylvain Chevallier

Joe Wiart

Isabelle Bloch

Time-frequency Selection in Two Bipolar Channels for Improving the Classification of Motor Imagery EEG*

Time and frequency information is essential to feature extraction in a motor imagery BCI, in particular for systems based on a few channels. In this paper, we propose a novel time-frequency selection method based on a criterion called Time-frequency Discrimination Factor (T FDF) to extract discriminative event-related desynchronization (ERD) features for BCI data classification. Compared to existing methods, the proposed approach generates better classification performances (mean kappa coefficient = 0.62) on experimental data from the BCI competition IV dataset IIb, with only two bipolar channels.

I. INTRODUCTION

Brain-computer interface (BCI) systems capture subject's intentions by measuring his brain activities and then translate them into commands to build a direct communication between brain and computer. Usually, brain activities are measured through electroencephalogram (EEG) because of its low cost and high time resolution [START_REF] Wolpaw | Brain-computer interfaces for communication and control[END_REF]. One of the most typical brain activities that can be captured by EEG recording is the band power modulation of sensorimotor rhythms, i.e. µ and β rhythms, during motor imagery [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF]. Existing studies show that imaginary movements of different body parts can cause the attenuation of sensorimotor rhythms, termed as eventrelated desynchronization (ERD), at corresponding "active" cortex areas; meanwhile, an enhancement of sensorimotor rhythms called event-related synchronization (ERS) might be observed at other "idling" areas [START_REF] Pfurtscheller | Mu rhythm (de)synchronization and EEG single-trial classification of different motor imagery tasks[END_REF]. Thus, motor imagery can be identified by classifying ERD/ERS patterns, which gives birth to motor imagery BCI [START_REF] Wolpaw | Brain-computer interfaces for communication and control[END_REF]. However, the poor signalto-noise ratio (SNR) of raw EEG signals and the overlap of different rhythms (e.g. α and µ rhythms) make it difficult to extract discriminative features for BCI classification [START_REF] Naeem | Dimensionality reduction and channel selection of motor imagery electroencephalographic data[END_REF]. Although multi-channel recording with a large number of electrodes (e.g. 64 or 128) and spatial filters can improve the SNR and extract discriminative features from overlapping signals [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis[END_REF], this setting reduces the portability and practicability of BCI so that it represents a main drawback for final users.

To solve this problem, several researchers proposed to record motor imagery EEG data with only a few bipolar electrodes placed around task-relevant sensorimotor areas of eaux, France. joe.wiart at orange-ftgroup.com the brain [START_REF] Lou | Bipolar electrode selection for a motor imagery based brain-computer interface[END_REF], [START_REF] Leeb | Brain-computer communication: motivation, aim, and impact of exploring a virtual apartment[END_REF]. Bipolar recording is an EEG acquisition technique that measures potential differences between two electrodes. This can improve the SNR by eliminating common additive noise of both electrodes (for details, see [START_REF] Lou | Bipolar electrode selection for a motor imagery based brain-computer interface[END_REF]). The positions of bipolar electrodes can be optimized algorithmically or using prior knowledge on the spatial location of brain activity during motor imagery (e.g. imagination of right hand movement elicits ERD around left sensorimotor areas). This allows reducing the number of electrodes while improving classification performances.

Beside spatial information, time and frequency information is also very important for classification, because motor imagery elicits ERD/ERS in specific bands within specific time segments instead of over all frequency bands in the whole trial. In general, µ rhythm is in 8-12Hz band and β rhythm in 18-25Hz band, but these rhythms can vary with subjects. Thus, a number of algorithms were proposed to find subject-specific frequency bands for extracting discriminative features [START_REF] Ang | Robust filter bank common spatial pattern (RFBCSP) in motor-imagery-based braincomputer interface[END_REF], [START_REF] Chin | Multi-class filter bank common spatial pattern for four-class motor imagery BCI[END_REF], [START_REF] Lemm | Spatio-spectral filters for improving the classification of single trial EEG[END_REF], [START_REF] Novi | Sub-band common spatial pattern (SBCSP) for brain-computer interface[END_REF], [START_REF] Yamawaki | An enhanced timefrequency-spatial approach for motor imagery classification[END_REF]. Most of them are based on extensions of the common spatial pattern (CSP) algorithm, such as filter bank CSP (FBCSP) [START_REF] Ang | Robust filter bank common spatial pattern (RFBCSP) in motor-imagery-based braincomputer interface[END_REF], [START_REF] Chin | Multi-class filter bank common spatial pattern for four-class motor imagery BCI[END_REF], subband CSP (SBCSP) [START_REF] Novi | Sub-band common spatial pattern (SBCSP) for brain-computer interface[END_REF] and common spatio-spectral pattern (CSSP) [START_REF] Lemm | Spatio-spectral filters for improving the classification of single trial EEG[END_REF], which combine flexible frequency filtering techniques with CSP to improve BCI performances. Among those methods, FBCSP seems to be the most effective one, because it yields the best BCI performances on multiple BCI competition data [START_REF] Ang | Robust filter bank common spatial pattern (RFBCSP) in motor-imagery-based braincomputer interface[END_REF], [START_REF] Chin | Multi-class filter bank common spatial pattern for four-class motor imagery BCI[END_REF]. Moreover, the robust FBCSP (RFBCSP), which includes robust Minimum Covariance Determinant (MCD) estimator in the algorithm, can be applied to motor imagery BCI with only three bipolar channels. However, RFBCSP involves mutual information (MI)-based feature selection and MCD estimator, both of which needing tedious iterative steps that greatly increase its complexity. Moreover, this method used a fixed time segment (i.e. 0.5-2.5s after cue on-set) of the EEG data for all subjects in the training process, which ignores the potential individual difference in optimal time segment selection.Although a recent work on FBCSP has involved optimal temporal selection process [START_REF] Ang | Mutual information-based selection of optimal spatial-temporal patterns for single-trial EEGbased BCIs[END_REF], it selected the time segment from only four different options (0.5-2.5s, 1.0-3.0s, 1.5-3.5s and 0.5-3.5s) and did not led to significant improvement on bipolar channel data (BCI competition IV IIb) compared to previous results in [START_REF] Ang | Robust filter bank common spatial pattern (RFBCSP) in motor-imagery-based braincomputer interface[END_REF].

In this paper, we propose a time-frequency selection method for extracting subject-specific ERD features for further improving motor imagery EEG classification with only a few bipolar electrodes. This method is based on a proposed criterion called Time-frequency Discrimination Factor (TFDF) to extract discriminative ERD features for the classification. Compared with the methods mentioned above, the proposed method is different in the following ways: 1) this method is neither based on CSP algorithm nor combined with other complex algorithms (e.g. MI-based algorithms or MCD estimator); 2) the subject-specific time-frequency area is selected from more than one thousand possible areas by a novel time-frequency selection criteria (TFDF), which can further improve the classification performances on bipolar recording BCI data; 3) this method can be applied on data recorded with only two bipolar channels for binary-class problems (e.g. a discrimination between imaginations of left and right hand movements).

II. EXPERIMENTAL DATA

We used data of the BCI competition IV dataset IIb [START_REF] Leeb | Brain-computer communication: motivation, aim, and impact of exploring a virtual apartment[END_REF] in this work. The dataset consists of two classes (left vs. right hand) cue-driven motor imagery BCI data from 9 subjects. The EEG data are recorded using 3 bipolar channels, i.e. at positions C3, Cz and C4. For each subject, 5 sessions are provided, including 3 training sessions and 2 evaluation sessions. The first two training sessions consist of 240 single trials (120 trials per session) without visual feedback. Each trial starts with a fixation cross and a short acoustic warning tone. Later a visual cue is given to guide the subject to execute the corresponding imagination of hand movement over a period of 4 seconds. The last training session (160 trials) and the two evaluation sessions (160 trials per session) are recorded with visual feedback from 0.5 to 4.5s after the cue on-set (for details, see [START_REF] Leeb | Brain-computer communication: motivation, aim, and impact of exploring a virtual apartment[END_REF]).

Many studies have shown that the electrode positions approximately overlaying the hand representation areas i.e. C3 (right hand) and C4 (left hand), usually generate the best results of single-trial discrimination between imaginations of left and right hand movements in non-CSP based classification [START_REF] Pfurtscheller | Mu rhythm (de)synchronization and EEG single-trial classification of different motor imagery tasks[END_REF], [START_REF] Pfurtscheller | EEGbased discrimination between imagination of right and left hand movement[END_REF], while Cz located at the functional area related to foot movement mainly contributes to the best discrimination between tasks involving imagination of foot movement [START_REF] Pfurtscheller | Mu rhythm (de)synchronization and EEG single-trial classification of different motor imagery tasks[END_REF]. Thus, we only used C3 and C4 in this study for improving the left v.s. right hand classification by using the proposed method.

III. TIME-FREQUENCY ANALYSIS AND SELECTION

A. Visualization of ERD/ERS in the Time-frequency Domain

The ERD/ERS patterns are usually expressed as percentage power decrease (ERD) or power increase (ERS) referring to the 1-s interval before the warning tone (for details, see [START_REF] Pfurtscheller | Mu rhythm (de)synchronization and EEG single-trial classification of different motor imagery tasks[END_REF]). The time-frequency maps of ERD/ERS for both left (L) and right (R) hands in the bipolar channels C3 and C4 were generated by the Biosig Toolbox using overlapping 2Hz bands (step = 1Hz) in the frequency range between 6 and 32Hz [START_REF] Schlogl | Biosig: a free and open source software library for BCI research[END_REF]. The obtained time-frequency maps are used to show the significant ERD/ERS (t-percentile bootstrap, α = 0.05) during the imaginations of right and left hand movements in C3 and C4, and to validate the necessity and effectiveness of the proposed time-frequency selection method in feature extraction (see Section IV).

B. Subject-specific Time-frequency Selection

Firstly, the signals in both channels are divided into EEG segments, in a series of overlapping time-frequency areas

{(ω m , τ n )} M×N with different frequency bands ω m = [ f m , f m + F -1], m = 1, ..., M, f m+1 = f m + F s (F
is the bandwidth, F s is the frequency step) and time intervals τ n = [t n ,t n + T -1], n = 1, ..., N,t m+1 = t m + T s (T is the interval width, T s is the time step).

As ERD features for the BCI classification are often measured by the logarithm of variance of band-pass filtered EEG in a specific time interval, so called logarithmic band power (BP) estimator [START_REF] Vidaurre | Time domain parameters as a feature for EEG-based brain-computer interfaces[END_REF], the overall BP of EEG segment (K trials) in a time-frequency area (ω m , τ n ) for each condition and each channel is estimated as follows.

Step 1: The variance of EEG segment in the time domain for each trial i, i = 1, ..., K for each condition and each channel is computed as:

v(i) = 1 T -1 t n +T -1 ∑ j=t n (x i j -xi ) 2 (1)
where x i j the j-th sample in the time interval τ n of the i-th trial of the ω m -bandpass filtered EEG data, and xi is the mean value over all samples of filtered EEG in the time interval τ n of the i-th trial.

Step 2: According to the existing literature, the overall BP for each class (L, R) and each channel (C3, C4) is calculated as the logarithm of median or mean of data variances over trials. Here, we use medians instead of mean values in the calculation, because outliers in the experimental data we used were reported, e.g. in [START_REF] Ang | Robust filter bank common spatial pattern (RFBCSP) in motor-imagery-based braincomputer interface[END_REF].

BP L C3 = log med v L C3 (i) (2)
BP L C4 = log med v L C4 (i) (3)
BP R C3 = log med v R C3 (i) (4)
BP R C4 = log med v R C4 (i) (5)
where med(•), denotes the median of the data over all trails, i = 1, ..., K, for each condition and each channel. Thus, the pattern difference (PD) between two conditions (left vs. right hand) in time-frequency area (ω m , τ n ) at each channel is expressed as:

PD C3 = BP L C3 -BP R C3 ( 6
)
PD C4 = BP L C4 -BP R C4 (7)
Imaginations of left and right hand movement usually elicited contrary contralateral dominance of ERD at C3 and C4 [START_REF] Pfurtscheller | Mu rhythm (de)synchronization and EEG single-trial classification of different motor imagery tasks[END_REF], [START_REF] Pfurtscheller | EEGbased discrimination between imagination of right and left hand movement[END_REF]. These task-spatial-related discriminative modulations can be measured by discriminative force Fd(ω m , τ n ) to estimate this positive contribution in a time-frequency area (ω m , τ n ) for the classification:

Fd(ω m , τ n ) = |PD C3 -PD C4 | = log med v L C3 (i) med v R C3 (i) / med v L C4 (i) med v R C4 (i) (8)
A large Fd(ω m , τ n ) indicates that large discriminative modulations occur in the time-frequency area (ω m , τ n ).

On the other hand, subjects are looking at the screen during both motor imagery tasks, which can generate some nonmotor-imagery-related (i.e. visual related) common modulations at C3 and C4 [START_REF] Blankertz | Optimizing spatial filters for robust EEG single-trial analysis[END_REF]. Meanwhile, neural activities at C3 and C4 may also affect each other due to volume conduction effects. These possible negative effects in the time-frequency area (ω m , τ n ) on the classification are estimated using a blurring force Fb(ω m , τ n ):

Fb(ω m , τ n ) = |PD C3 + PD C4 | = log med v L C3 (i) med v R C3 (i) • med v L C4 (i) med v R C4 (i) (9)
A small Fb(ω m , τ n ) indicates that small common modulations happen in the time-frequency area (ω m , τ n ).

We call the difference between Fd(ω m , τ n ) and Fb(ω m , τ n ) Time-frequency Discrimination Factor, T FDF(ω m , τ n ), and use it to evaluate the overall contribution of the data in the time-frequency area (ω m , τ n ) from electrodes C3 and C4 for the binary-class discrimination:

T FDF(ω m , τ n ) = Fd(ω m , τ n ) -Fb(ω m , τ n ) ( 10
)
An ideal time-frequency area for classification should have large discriminative modulations and small common modulations, so that the optimal time-frequency area (ω , τ opt ) is estimated by seeking the maximum value of T FDF(ω m , τ n ) among all time-frequency areas.

C. Data Elaboration, Feature Extraction and Classification

In this study, the sessions used for training are selected based on the suggestions of winners in BCI competition IV [START_REF] Blankertz | BCI competition IV[END_REF] in order to compare our results with theirs. For each bipolar channel, 5th order butterworth filters are applied to compute 19 successive 4Hz-wide frequency bands of signals: 8-12Hz, 9-13Hz, 10-14Hz, ..., 26-30Hz, and 15 successive 8Hz-wide frequency bands of signals: 8-16Hz, 9-17Hz, 10-18Hz, ..., 22-30Hz (M = 34). Then 39 overlapping time segments in each frequency band were obtained through 2s, 2.5s and 3s-wide sliding windows (13 segments for each sliding window), respectively, with 0.2s-step moving from 0.5s after the cue on-set (N = 39). Those parameters are set based on the experience from competitors reported in BCI competition IV [START_REF] Blankertz | BCI competition IV[END_REF]. Therefore, there are 34 × 39 = 1326 time-frequency areas for subject-specific selection. The BP values are extracted from two channels of EEG in the selected time-frequency area (ω opt , τ opt ) as ERD features. Then, Fisher's linear discriminant analysis (LDA) was employed with these features for the classification as most BCI classification works did [START_REF] Lotte | A review of classification algorithms for EEG-based brain-computer interfaces[END_REF].

IV. EXPERIMENTAL RESULTS AND DISCUSSION

Fig. 1 shows an example of time-frequency maps displaying significant ERD (red) and ERS (blue) for a typical subject (Subject 6) in BCI competition IV dataset IIb [START_REF] Leeb | Brain-computer communication: motivation, aim, and impact of exploring a virtual apartment[END_REF]. From Fig. 1, we can see that 1) imagination of left hand movement elicits C4 dominant µ band ERD patterns while that of right hand movement generates C3 dominant µ band ERD patterns around 12Hz approximately during 0.5-5s after cue on-set; 2) a background noise (generated by non-motorimagery-related activities) is widely observed in the timefrequency map (e.g. in 15-30Hz, 1-4s) so that information from those areas may deteriorate the classification results. Based on these observations, we can conclude that 1) timefrequency selection is necessary for extracting discriminative ERD patterns for a good classification performance; 2) common noise cannot be completely removed in all timefrequency areas through bipolar recoding technique. Thus, a subject-specific time-frequency selection method is required to capture the ERD features in the two bipolar channels that contain discriminative information. Table I lists the optimal frequency band ω opt and time interval τ opt of the time-frequency areas selected from the training data by the proposed method for each subject. For example, the optimal frequency band of Subject 6 is 10-14Hz and the optimal time segment is 0.7-2.7s, which is slightly different from the general µ band (8-12Hz) and the suggested time segment (0.5-2.5s) in the literature [START_REF] Ang | Robust filter bank common spatial pattern (RFBCSP) in motor-imagery-based braincomputer interface[END_REF]. From Fig. 1, we can see that this time-frequency area contains the discriminate ERD patterns of two classes: imaginations of left hand and right hand movements elicited contrary contralateral ERD dominance in the selected areas, which indicates the effectiveness of the proposed method in timefrequency selection for discrimination. The optimal time segments and frequency bands for capturing discriminate ERD patterns vary with different subjects (see Table I), so that it is necessary to use subject-specific time-frequency area in the classification.

In the classification, the EEG signals in the subjectspecific time-frequency area (ω opt , τ opt ) from all training trials are used for training, then the ω opt -bandpass filtered EEGs segments with the same time length of τ opt (i.e. T ) are obtained from each entire single-trial of testing data via a

s) 2.9-4.9 2.7-5.2 0.5-2.5 0.9-2.9 0.7-3.2 0.7-2.7 1.1-3.1 0.9-2.9 0.9-2.9 0.2s step sliding window to generate continuous classification outputs. According to the BCI competition requirement, the classification performance is measured by kappa coefficient κ value [START_REF] Dornhege | Toward brain-computer interfacing[END_REF]. Table II provides the mean performance over the two independent evaluation sessions obtained by our method, and comparisons with the three best submissions of the BCI competition IV on this dataset [START_REF] Blankertz | BCI competition IV[END_REF]. Our method generates the best mean kappa value ( κ = 0.62) among all methods on these data. To further analyze the results, the paired-sample t-test was employed to reveal the statistically significant difference between the performances of different methods. Although improvements of kappa values compared 1 winner ( κ = 0.60, p = 0.12), who used RFBCSP [8] and 2 nd winner( κ = 0.58, p = 0.18), who applied common spatial subspace decomposition (CSSD) on multiply frequency bands [START_REF] Blankertz | BCI competition IV[END_REF], are not statistically significant, our method outperforms the 1 st winner in 6 out of 9 subjects (except subjects 4, 7 and 8), and the 2 nd winner in 6 out of 9 subjects (except subjects 4, 7 and 9). For this dataset, our results are always better than the 3 rd placed one ( κ = 0.46, p < 0.01), who employed classical CSP with spectrally filtered neural time series prediction preprocessing (NTSPP) signals [START_REF] Blankertz | BCI competition IV[END_REF]. As all of their works have used all three bipolar channels (C3, Cz and C4) provided by the dataset, our method not only generated better mean performances but also used less channels than them, which shows the interest of our method in feature extraction and electrode reduction.

V. CONCLUSION

Previous state of the art has highlighted the need for user-specific parametrization in the preprocessing steps of motor imagery BCI, especially in term of frequency band adaptation. This contribution, through the proposition of a novel algorithm and a new criterion, emphasizes that the calibration in the temporal domain is also relevant. The results show that our approach yields better average performance using only two bipolar channels and a simple LDA classifier compared to those obtained by existing methods on the BCI competition dataset IIb. However, this conclusion has to be confirmed on more datasets; and a deeper study of the differences with respect to existing methods should be conducted. In the future, we will also extend this method to multi-class problems.

Fig. 1 .

 1 Fig. 1. Time-frequency maps displaying significant ERD (red) and ERS (blue) for Subject 6 (a typical example) in BCI competition IV IIb. The areas in the rectangles are the time-frequency areas selected by the proposed method. Imaginations of left and right hand movement elicited contrary contralateral dominant ERD in the selected areas, which indicates the effectiveness of the proposed method in time-frequency selection for discrimination.

TABLE I FREQUENCY

 I BAND ω opt AND TIME INTERVAL τ opt SELECTED BY THE T FDF VALUES FOR ALL THE SUBJECTS.

					Subjects
	1	2	3	4	5	6	7	8	9
	ωopt (Hz) 10-14 11-15	9-13	8-12	26-30 10-14 11-15 10-14 22-26
	τopt

TABLE II COMPARISON

 II OF THE CLASSIFICATION PERFORMANCES BETWEEN THE PROPOSED T FDF METHOD AND OTHER METHODS. st ) 0.40 0.21 0.22 0.95 0.86 0.61 0.56 0.85 0.74 0.60 CSSD (2 nd ) 0.42 0.21 0.14 0.94 0.71 0.62 0.61 0.84 0.78 0.58 CSP+NTSPP(3 rd ) 0.19 0.12 0.12 0.77 0.57 0.49 0.38 0.85 0.61 0.46

		Subjects
		1 2 3 4 5 6 7 8 9 Mean
	T FDF	0.44 0.24 0.25 0.93 0.86 0.70 0.55 0.85 0.75 0.62
	RFBCSP(1

*This work was partially supported by grants from China Scholarship

Council and Orange Labs.


Estelle Duveau

Simon Courtemanche

Lionel Reveret

Edmond Boyer

Cage-based Motion Recovery using Manifold Learning

We present a flexible model-based approach for the recovery of parameterized motion from a sequence of 3D meshes without temporal coherence. Unlike previous model-based approaches using skeletons, we embed the deformation of a reference mesh template within a low polygonal representation of the mesh, namely the cage, using Green Coordinates. The advantage is a less constrained model that more robustly adapts to noisy observations while still providing structured motion information, as required by several applications. The cage is parameterized with a set of 3D features dedicated to the description of human morphology. This allows to formalize a novel representation of 3D meshed and articulated characters, the Oriented Quads Rigging (OQR). To regularize the tracking, the OQR space is subsequently constrained to plausible poses using manifold learning. Results are shown for sequences of meshes, with and without temporal coherence, obtained from multiple view videos preprocessed by visual hull. Motion recovery applications are illustrated with a motion transfer encoding and the extraction of trajectories of anatomical joints. Validation is performed on the Hu-manEva II database.

Introduction

Motion recovery using multiple view environments has known important developments in Computer Vision and Graphics. Applications to biomechanics and character animation require a reliable and interpretable measurement of motion. In that respect, the recovery of articulated motion appears as a relevant goal since the underlying skeleton provides a deterministic interpretation of motion. As a preprocess, it is now common to calibrate multiple views to obtain streams of 3D meshes from visual inputs such as silhouettes typically using space carving or visual hull [START_REF] Corazza | Automatic generation of a subject-specific model for accurate markerless motion capture and biomechanical applications[END_REF]. The resulting information from visual hull is not temporally structured yet nor related to interpretable parameters in terms of motion. In contrast to existing methods for mesh tracking that use the coupling of an articulated skeleton and skinning to deform a known template, e.g. [START_REF] Vlasic | Articulated mesh animation from multi-view silhouettes[END_REF][START_REF] Urtasun | 3d people tracking with gaussian process dynamical models[END_REF], we propose here to relax the strong parametric constraint of the skeleton by making use of barycentric coordinates. In barycentric coordinates, a low polygonal version of the reference template, the cage, controls the deformation of the mesh by interpolation. We use Green Coordinates for their known properties of good conformal behavior [START_REF] Lipman | Green coordinates[END_REF]. The output is thus a sequence of temporally coherent coordinates of the vertices of the cage, controlling the deformation of a mesh template and constrained to lie on a pose manifold. To fulfill the goal of motion recovery, the cage has been carefully designed to follow a morphologically standardized representation of the body shape. It can thus be typically converted to an articulated skeleton representation.

For most methods, the accuracy of articulated tracking may still be challenged by lighting conditions, poor camera positioning or occlusions. One of the key factor is therefore the ability to bring an efficient regularization to the tracking problem. Latent spaces have been successfully used to constrain the manifold of plausible poses for human body tracking from monocular video or pose editing in animation. In this paper we investigate a similar strategy but in the novel context of motion recovery as a postprocess of 3D reconstruction of meshes from multiple view silhouettes. The direct link between the morphological cage and an articulated skeleton allows the use of available databases of motion capture data for the automatic learning of the manifold.

Our main contribution is thus the derivation of a mesh tracking algorithm, based on a Green Coordinates encoding, and regularized by a manifold learned whether by hand labeling of poses or automatically from motion capture data. We build on the formalism introduced by Prisacariu and Reid [14][15] for a rigorous derivation of the gradients of the object measurement with respect to the space of latent variables. We also provide the exact formulation of the gradients. We demonstrate our approach on different human performances : on temporally coherent clean meshes for consistency test [START_REF] Vlasic | Articulated mesh animation from multi-view silhouettes[END_REF], on the HumanEva II dataset [START_REF] Sigal | Humaneva: Synchronized video and motion capture dataset and baseline algorithm for evaluation of articulated human motion[END_REF] for objective evaluation and on non temporally coherent meshes obtained from a personal set-up adapted to visual hull for biometrics test.

Related Work

Motion recovery from a sequence of meshes

For mesh tracking, feature-level methods such as the patch-based approach of [2] offer a maximum flexibility. This approach is robust to large changes in the topology of the analyzed scene as no assumption is made on the internal structure but it does not fulfill the objective of motion recovery as no interpretable parametrization can be directly derived from this formulation. Following the parametrization encoded by our normalized cage, we deliver a more structured output. Furthermore, in case of very noisy input meshes, an approach that is too local risks providing unreliable results with non recoverable drift. This effect appears clearly on the HumanEva dataset [START_REF] Sigal | Humaneva: Synchronized video and motion capture dataset and baseline algorithm for evaluation of articulated human motion[END_REF], where the camera locations are not optimal for a preprocess step using visual hull reconstruction, leading to noisy input meshes.

Taking an opposite direction to feature-level approaches, the skeleton-based methods impose a strong prior on the parameter space. This approach has shown successful results such as in [START_REF] Vlasic | Articulated mesh animation from multi-view silhouettes[END_REF] and [8]. In both cases, a skeleton of a human or animal character has to be aligned on a reference template mesh. A skinning algorithm deforms the mesh, which is matched against 2D silhouettes. We propose here an alternative method for the deformation phase using Green Coordinates [START_REF] Lipman | Green coordinates[END_REF] and a 'cage' proxy to control the shape. The skinning algorithm is known to be sensitive to a careful tuning of vertices weights, where the automatic phase usually needs for manual adjustment. The vertices weights automatically provided by the Green Coordinates method show superior stability and do not require further tuning.

Tracking mesh data using a cage for barycentriccoordinates has first been proposed by [START_REF] Savoye | Cage-based Tracking for Performance Animation[END_REF] and recently improved by [START_REF] Thiery | Cager: From 3d performance capture to cage-based representation[END_REF]. Unlike these works, we make no assumption on the temporal coherence of the input 3D data. Another major difference is in the regularization strategy. In [START_REF] Savoye | Cage-based Tracking for Performance Animation[END_REF], the cage deformation was constrained by a Laplacian smoothing. It has been reported that it did not prevent the cage from collapsing when tracking long sequences. This is our motivation for investigating a more robust regularization strategy using manifold learning.

Manifold learning for human motion tracking

The success of manifold learning for robust tracking of human motion from monocular videos has been typically proved with works such as [START_REF] Urtasun | 3d people tracking with gaussian process dynamical models[END_REF] which uses GPLVM. For animation, an application of manifold learning has also been shown for 3D pose editing [9]. Other methods for animation editing using low dimensional subspace embedding exist. Among recent works, [START_REF] Cashman | A continuous, editable representation for deforming mesh sequences with separate signals for time, pose and shape[END_REF] uses different encodings for shape, pose and time. However, it requires an animation sequence as input whereas our learning samples are few and sparse. In our case, we therefore use a similar to [9] manifold learn-ing approach, but in the novel context of cage-based animation from a sequence of meshes. We build on [START_REF] Prisacariu | Nonlinear shape manifolds as shape priors in level set segmentation and tracking[END_REF], who uses GPLVM to segment 2D images, by using a non-linear manifold as a deformation prior to track a cage-based animation. In particular, closed-form solutions of gradients are provided for a straightforward implementation of a simple optimization based on gradient descent.

The rest of the paper is organized as follows : in section 3, we describe how our template is defined and deformed. In section 4, we explain the learning procedure to obtain a latent space of coherent motion. In section 5, we detail latent space-based tracking. Section 6 shows results obtained with our method and evaluation.

Deformation Model

For our model-based approach, we need a shape model and a deformation model. Our shape model is a 3D mesh M. To deform this mesh, we use a cage-based approach as it has better properties in terms of shape preservation and smoothness than skinning. In addition, it makes the motion parameters independent on the actual model as the same cage is used for every specific human model. Our template surface is defined as a mesh embedded in a coarse bounding mesh. The template mesh is called the model, the coarse bounding mesh is called the cage. The barycentric coordinates used to generate the model from the cage are given by the Green coordinates [START_REF] Lipman | Green coordinates[END_REF]. As mentioned in [START_REF] Lipman | Green coordinates[END_REF], unlike other barycentric coordinates, Green coordinates are empirically shown to be quasi-conformal in 3D. There are no requirements on the model : the only information used is the position of its vertices. Green Coordinates behave properly if there are no self-intersections in the cage when computing the coordinates. However, once the coordinates are computed, self-intersecting deformations of the cage do not create artefacts. Model vertices are expressed as a linear combination of cage vertices V = (v 1 . . . v n ). Therefore a point P of the model is deformed according to :

P → n i=1 φ i (P )v i + m j=1 ψ j (P )n(t j ) (1)
where (t 1 . . . t m ) are the faces of the cage, n(t j ) is the normal of face t j and {φ i } i=1...n , {ψ j } j=1...m are the Green coordinates.

Being a coarse approximation of the model, the cage is easy to deform. However, it lacks an underlying structure to intuitively generate and deform it in a volume and shapepreserving way. Works such as [START_REF] Ben-Chen | Variational harmonic maps for space deformation[END_REF] allow the user to interact with the cage using a small number of position and orientation constraints. However, how to generate the cage remains an issue. Inspired by skeletal joints, in order to provide a standard way of building and deforming the cage, we define oriented 3D quads located at main joints. These quads Each oriented quad has a constant size which is not modified during the motion. The cage is then built procedurally from the vertices defined by the quads. Deforming the model then consists in setting the translation and orientation of each quad, which procedurally defines the vertices of the cage. The topology of the cage is computed once by defining relations between quads. The geometry of the cage is updated every time the OQR is modified with the positions of the vertices of the quads. For the sake of brevity, details could not be included in the present paper, which focuses on motion recovery. Finally, this cage defines the vertices of the model by linear combination of the vertices and normals of the cage (see Equation ( 1)).

As a result, the number of degrees of freedom of the deformation is greatly reduced. From 3 times the number of vertices of the model, we reduce it to 6 times the number of quads. In the example on Figure 1, it is reduced from 3 × 10002 to 6 × 21. This approach can be interpreted as an intermediate representation of an articulated character, tween a hierarchical skeleton bind to a mesh using skinning and a standard non-hierarchical mesh. For the remainder of the paper, we refer to this representation as Oriented Quads Rigging (OQR).

Learning a manifold of deformations

Thanks to the constant size of the quads, some unrealistic deformations such as shrinking the circumference of an arm are avoided. However, nothing preserves the length of the limbs or non physiological configuration. The learning of a non-linear manifold of the possible positions and orientations of the quads allows to enforce such constraints. As in [START_REF] Prisacariu | Nonlinear shape manifolds as shape priors in level set segmentation and tracking[END_REF], where Gaussian Process Latent Variable Models are used to segment 2D images, we use GPLVM to learn a manifold of deformations. A Gaussian Process Latent Variable Model (GPLVM) is a low dimensional representation of the original data called the latent space. Each data point y i of dimension d can be represented by a latent point x i of dimension q, where q is smaller than d. Conversely, every Figure 2: Example of 2D-latent space and corresponding deformations. The gray level corresponds to the variance at a given point in latent space. Blue crosses are the latent coordinates of the learned data-set. Red points and their corresponding reconstructions show that the latent space is suitable for interpolation between learned poses. point x in the latent space projects to a data point through a Gaussian distribution defined by its mean µ and its variance σ 2 . As derived in [START_REF] Lawrence | Probabilistic non-linear principal component analysis with gaussian process latent variable models[END_REF], for a data set of size n, we have :

P (Y |X, θ) = n i=1 N (y i |0, K(θ)) (2)
where Y = [y 1 . . . y n ], X = [x 1 . . . x n ] and K i,j = k(x i , x j ) with k being the covariance function defined by the hyper-parameters θ. For instance, if k is a radial basis function, as in our case, it is written :

k(x i , x j ) = αexp(- γ 2 (x i -x j ) T (x i -x j )) (3)
and the hyper-parameters are α and γ. To learn the GPLVM, Equation (2) is maximized with respect to both X and θ.

In our case, each data point y i describes one deformation by the positions and orientations of the quads in the OQR. We align all the shapes with respect to the position and orientation of the quad at the pelvis. We then learn the latent variables and hyper-parameters by maximizing Equation (2). An example of a 2D-latent space can be seen on Figure 2. We impose back-constraints so that points close in data space stay close in latent space, which is not guaranteed in the usual formulation of GPLVM [12].

Tracking

To recover the motion, we must fit the model to the current representation of the shape, namely the target, which is, in our case, a 3D mesh. To do so, we derive a differentiable energy that measures how well the model fits the target, the model being the template mesh controlled by the cage. We differentiate this energy in the latent space to recover the deformation of the model. By doing so, we stay in the learned manifold i.e. in the space of plausible shapes.

Distance between two meshes

We make no assumptions on the temporal coherence of the target meshes. To measure how well the model mesh fits the target mesh, we thus need to measure the distance between two different meshes having potentially two different topologies. As we will minimize this distance, we need a function which is differentiable with respect to the model mesh. To that end we use a distance field built such that the target mesh is its zero level-set. Distance fields have the advantages of being easy and fast to compute. They are differentiable everywhere and can be used to measure the distance between meshes of different topologies.

Other approaches exist to deal with distances between non-corresponding meshes. One can first estimate vertexto-vertex sparse correspondences, then use those pairs of vertices to compute the distance between the two meshes. However, computing correspondences is costly and can fail. Another approach is to directly use a mesh to mesh distance, such as the Hausdorff distance. However, such distances are often hard to compute and can not be differentiated because it implies finding optima. We thus use a distance field to avoid these issues [10]. Figure 3 shows examples of slices of distance fields that we obtain on noisy visual hulls.

Once the distance field is computed, we can define the distance function that has to be minimized. For a given mesh M, the distance to the target T is computed as the sum of the distances of each vertex of M to T . The distance of a vertex to T is the value of the distance field at the voxel to which the vertex belongs. Thus the distance function takes the following form :

d DF (M, T ) = Pi∈M d DF (P i , T ) (4)
where d DF (P i , T ) is the value of the distance field generated from T at position P i .

Shape Optimization

To recover the deformation of the model, we could minimize the distance between the model and the target in the unconstrained highly-dimensional OQR space. To limit the evolution of the model to plausible deformations, we differentiate this distance in the latent space. As a result, we stay in the low-dimensional learned manifold i.e. in the space of plausible poses.

We define an energy function to minimize with two terms. The first one measures how close the model and the target are. This is given by their distance as defined in Section 5.1. The second term takes advantage of the fact that at each point in the latent space, we can compute the variance which indicates how likely this latent point is to generate a valid data point. We use this variance as a regularization term to ensure that the shape generated from the latent point is likely to exist. We therefore try to find :

x = arg min x (d DF (M(x), T ) -log(σ 2 (x))) (5)
where x is the latent point, M(x) is the model generated from the latent point x and σ 2 is the variance at point x.

We do this by differentiating this energy with respect to the latent variables. This is done by using the chain rule :

d(d DF (M, T )) dx = d(d DF (M, T )) dM dM dV dV dy dy dx (6)
where V is defined in Section 3 as the vertices of the cage and y is the deformed shape generated from x i.e. the set of positions and orientations of the oriented quads. Each vertex in V is connected rigidly to one and only one oriented quad. dV dy follows trivially. y is the mean of the Gaussian distribution by which the latent point x projects to a data point :

y = Y T K -1 K x ( 7
)
where K is the covariance in-between all the learned latent points X = [x 1 . . . x n ] and K x is the covariance between x and all the learned latent points. Differentiating Equation [START_REF] Gall | Optimization and filtering for human motion capture[END_REF] gives :

dy dx = Y T K -1 dK x dx (8)
Each element of K x is the covariance function k applied on x and one learned latent point. As a result, dKx dx follows from Equation (3).

As for the regularization term, x projects to a data point with a mean as expressed in Equation ( 7) and a variance σ 2 such that :

σ 2 = k(x, x) -K T x K -1 K x (9)
As a result :

dσ 2 dx = dk(x, x) dx -2K T x K -1 dK x dx (10)
Thanks to Equations ( 6) and ( 10), a simple first order method such as gradient descent can be used to minimize Equation ( 5) and recover the deformation of the model.

Pose Optimization

So far, we have captured the deformation of the shape as encoded in the latent space. As explained in Section 4, this shape is defined up to the global transformation Γ of the pelvis. As a result, we also need to recover the global translation and orientation of the model.

To do so, we minimize the energy without the regularization term by differentiating Equation ( 4) with respect to the pose parameters :

d(d DF (M, T )) dΓ = d(d DF (M, T )) dM dM dΓ ( 11
)
Pose and shape optimization can be done either jointly or iteratively until convergence. Joint optimization is more correct but converges more slowly as the search space is larger. Empirically, iterating between pose and shape optimization works faster and better, especially when tracking a sequence. Indeed, the rigid transformation and the deformation between two frames are small enough that optimizing first for the pose does not lead to a pose from which finding the good shape is impossible, and vice-versa.

Experiments and Discussion

We tested our method on three test sets. We first tested our approach on a temporally coherent mesh of high quality. We then tackled more challenging data by tracking a sequence of noisy visual hulls, without temporal coherence, learning the motion manifold automatically from the CMU Motion Capture Database (http://mocap.cs.cmu.edu/). This was first evaluated on the HumanEva II dataset [START_REF] Sigal | Humaneva: Synchronized video and motion capture dataset and baseline algorithm for evaluation of articulated human motion[END_REF]. We finally used it for gait analysis. All experiments have been done with the deformation model of Figure 1. It consists in a rigging of 21 oriented quads, which generates a cage of 86 vertices and 166 faces.

Evaluation on a temporally coherent sequence

We first tested our deformation optimization without the pose recovery. To do so, we use a bouncing sequence. The sequence of target meshes is made of a temporally coherent animated mesh of 10000 vertices and 20000 faces. We define the template mesh as one of these meshes simplified to 4000 vertices and 8000 faces. The training was done by manually locating the oriented quads on 16 frames of the sequence. We learn a 10-D manifold from these samples. Figure 4 shows the progression of the gradient descent towards an optimal deformation. We can notice that, since we are optimizing in the latent space, each of the sample deformations is a valid deformation. Figure 5a shows the advantage of tracking in the manifold. As there is no structure on the cage, optimizing with respect to the oriented quads creates a shrinking of the cage as the quads are attracted to the closest patch of surface without any influence from the position and orientation of the other oriented quads. We can see on Figure 5 the influence of the regularization term of the energy. Without it, the model mesh is attracted to the shape that is the closest to the target mesh and can be generated from the latent space without any regard for the likelihood of this shape. Using the variance as a regularization term forces the optimal latent point to be more likely to generate a valid deformation. As explained in Section 5.3, we also recover the pose. In practice, we noticed that pose recovery and deformation re- covery done iteratively works better than joint optimization. Figure 6 shows results on the bouncing sequence.

Automatic learning on CMU database and validation on the HumanEva II dataset

Automatic learning : Instead of creating the training data by manually fitting the OQR to meshes, we also propose a method that uses existing databases of movements. From walking sequences in the CMU Motion Capture Database (http://mocap.cs.cmu.edu/), we automatically generate 120 training samples by rigidly mapping the quads of the OQR to the skeleton. It has to be done only once for a given frame of one of the subject's performance. After that, more samples can be automatically added to the training set using any available motion capture data samples for this subject. This training data is then used to build a 2D-manifold.

HumanEva II : We use the walk sequence of subject S4 in the HumanEva II dataset [START_REF] Sigal | Humaneva: Synchronized video and motion capture dataset and baseline algorithm for evaluation of articulated human motion[END_REF] for a quantitative evaluation. This data consists of 4 synchronized color cameras.

Using the background subtraction software and the mixture of Gaussian distribution model provided with the dataset, we extracted the silhouettes from the videos. From these silhouettes, we generated the visual hulls using the Exact Polyhedral Visual Hull method [6]. As the scene was captured with 4 cameras facing each other, the visual hulls are a very coarse approximation of the observed shape as can be seen on Figure 7a. The model used for cage-based deformation is the mesh model for subject S4 generated by Stefano Corazza and the Stanford BioMotion Group using [START_REF] Corazza | Automatic generation of a subject-specific model for accurate markerless motion capture and biomechanical applications[END_REF].

Local adjustment : One limitation of our method is its dependence on the training data. As with every manifoldbased method, we are limited to what the manifold can reconstruct. To account for this limitation, we can first optimize with respect to latent variables to get as close as possible to the target mesh. Once this has converged, a few more iterations can be done with respect to the oriented quads to locally adjust the cage by differentiating the energy function with respect to the oriented quads : An example of local adjustment is shown in Figure 8, where a 2D latent space does not capture enough variance to track the mesh properly. Table 1 shows the 3D error as specified in the HumanEva II dataset with and without the local adjustment for the first 150 frames. As our approach is not skeleton-based, the positions of the joints are computed by cage-based animation. On the first frame where motion capture data is available i.e. where the positions of the joints are known, the Green Coordinates are computed on these positions and then used to compute the new positions of the joints when the cage is deformed. With local adjustment, the mean error drops from 90mm to 38.6mm.

Figure 9 shows reprojections of the model tracked on the 7a) to reach the levels of tracking without local adjustment. sequence. We refer to the video for the whole sequence.

For the whole walking sequence, the error is 85.2mm (Table 1). However, when considering only the first 150 frames, where the visual hulls are of better quality (Figure 7b), the error drops to 38.6mm. The error for the whole sequence is about the same as the error for the first 150 frames without local adjustment. This can be explained by the non-discriminative quality of the target mesh. Indeed, the coarseness of the visual hull (Figure 7a) prevents any fine adaptation of the shape by local adjustment. Table 1 also reports the results found in [START_REF] Corazza | Markerless motion capture through visual hull, articulated icp and subject specific model generation[END_REF] and [START_REF] Gall | Optimization and filtering for human motion capture[END_REF]. The method of Corazza et al. [START_REF] Corazza | Markerless motion capture through visual hull, articulated icp and subject specific model generation[END_REF] also takes as input a series of visual hulls. The results are limited to the first 150 frames as reinitialisation is necessary after that due to the poor quality of the visual hulls. Thanks to the use of a manifold of deformations, even though the error increases after frame 150, we are still able to attain good results without re-initialisation.

To our knowledge, Gall et al. [START_REF] Gall | Optimization and filtering for human motion capture[END_REF] report the best results. However, they also use information not taken into consideration in this paper such as the appearance of the model.

Applications

Motion analysis: For gait analysis, walking and jogging sequences are selected in the CMU motion capture database Examples of tracking are shown on Figure 10. From this sequence, we can extract locomotion cycles. Figure 11 shows the height of the feet and the swinging motion of the arms during a cycle. Because of the noise from the treadmill that can be seen on Figure 10, the trajectories of the feet are slightly noisier that those of the hands. We refer to the video for the results on varying speeds of walks and jogs.

Motion transfer : Once we have captured our parameterized performance, we can put any model we want in the cage to obtain a new animation (motion transfer). We can edit the original shape or put a texture on it as there is no drift in vertices position. This process is fully automatic as the computation of the Green Coordinates is automatic and gives smooth deformations. On the contrary, skinning usually requires user intervention to correct skinning weights, which makes the process of changing the model more tedious. Figure 12 shows an example of the transfer of the tracking shown in Figure 6 to a different 3D model.

Conclusion

In this paper, we tackle the problem of motion recovery from a sequence of 3D meshes. To this end, we have proposed a method using a cage-based deformation constrained by manifold learning. The use of a manifold regularizes the tracking problem, leading to a more robust solution. We use oriented quads as an underlying structure for cage-based animation. This gives an intuitive way to generate and deform cages, the Oriented Quads Rigging (OQR). On top of that, the use of cage-based animation allows for the tracking of any model, not just models with an intrinsic skeletal structure. Experimental results demonstrate the validity of our approach. We have evaluated our method on temporally coherent data and tested its usability by generating the training data from an existing motion capture database. However, our current method has not considered the fact that classic GPLVM struggles with learning a good model for complex datasets and becomes intractable for large datasets. Using a stochastic approach as in [START_REF] Yao | Learning probabilistic non-linear latent variable models for tracking complex activities[END_REF] may overcome these issues and lead to a better tracking on complex sequences.

Figure 1 :

 1 Figure 1: The Oriented Quads Rigging (OQR): template model and example of deformation.

Figure 3 :

 3 Figure 3: Example of distance fields -model and slicing plane (left), corresponding slice of signed distance field (right) as well as another slice in the sagittal plane.

d

  (d DF (M,T )) dM is computed by finite differences.

dM dV is computed by differentiating Equation ( 1 )

 1 with respect to the vertices of the cage.

Figure 4 :Figure 5 :

 45 Figure 4: (a) Target model -(b) Initial pose -(c) Samples on the convergence trajectory.

Figure 6 :

 6 Figure 6: Optimisation of both pose and deformation. The tracking (gray) is overlaid over the target (yellow).

Figure 7 :Figure 8 :

 78 Figure 7: Examples of tracking on the HumanEva II dataset : (from left to right) visual hull, result and overlay of the model on the visual hull. (a) The location of the cameras can lead the visual hull to be a very coarse approximation of the observed shape. (b) When the visual hull quality increases, the tracking becomes more precise.

Figure 9 :

 9 Figure 9: Reprojection on camera 3 of the tracked model for frames 5, 45, 102 and 148 of the sequence for subject S4 in the HumanEva II dataset.

Figure 10 :

 10 Figure 10: Pose and shape recovery -The tracking (gray) is overlaid over the target (yellow). The training data is built automatically from the CMU motion capture database and the target meshes are independent visual hulls.

Figure 11 :

 11 Figure 11: Application : motion analysis -(a) Absolute height of the right and left ankles -(b) Position relative to the pelvis in the sagittal plane of the right and left wrists.

Figure 12 :

 12 Figure 12: Application : motion transfer from Figure 6.

Table 1 :

 1 Error for sequence of subject S4 of the HumanEva II dataset. The mean, median and standard deviation in mm is measured for the first 150 frames without (third row) and with (fourth row) the local adjustment. The error increases after frame 150 (bottom row) as the quality of the visual hulls decreases (Figure

	Method	Frames	Mean error (mm)	Median error (mm)	Standard deviation (mm)
	Corazza [5]	1-150	80.0		13.0
	Gall [7]	2-1258	32.01		4.53
	Our method -no local adjustment	1-150	90.0	87.3	24.1
	Our method	1-150	38.6	39.2	4.1
	Our method	1-300	85.2	55.8	55.8

Acknowledgments This work has been partially funded by grant ANR-2010-BLAN-0206-01 MORPHO of the French National Research Agency.