Phat Luong

Nguyen

Julien Mille

Dominique H Li

Donatello Conte

Nicolas Ragot

Efficient Dynamic Texture Classification with Probabilistic Motifs

We propose to tackle dynamic texture video classification as a pattern mining problem. In a nutshell, videos are represented by frequent sequences of representative patches. Firstly, we use a Gaussian Mixture Model to make the clustering of patches from training videos. Secondly, a soft assignment is used as an encoding method to construct sequences of probability vectors (p-sequences) representing sequences of spatio-temporal patches. Thirdly, for each class, we mine meaningful motifs appearing inside the training p-sequences by means of an adapted data mining approach. Finally, feature vectors are constructed from the mined motifs, using the probabilistic support, which quantifies the match between the p-sequences, of the video to be classified, and the key-motifs of the classes. Experimental results and analysis for dynamic texture classification on benchmark datasets (i.e. UCLA, Traffic) show the interest of the proposed method.

I. INTRODUCTION

Dynamic textures are repetitive spatio-temporal patterns characterized by non-rigid and complex motions [START_REF] Soatto | Dynamic textures[END_REF]. Extraction of spatio-temporal features and patterns may allow the characterization and classification of dynamic texture videos. In this paper, we are interested in creating a new representation of dynamic textures and measure its robustness by means of video classification, with a view to provide more insight on discriminant features than what can be expected from deep learning approaches. Let us first review existing texture video classification methods, i.e. hand-crafted features extraction and deep learning methods.

Early hand-crafted spatio-temporal feature-based methods depend on optical flow [START_REF] Nelson | Qualitative recognition of motion using temporal texture[END_REF]- [START_REF] Crivelli | Motion textures: modeling, classification, and segmentation using mixed-state markov random fields[END_REF]. For example, in the work of [START_REF] Crivelli | Motion textures: modeling, classification, and segmentation using mixed-state markov random fields[END_REF], a motion/no-motion map is built and combined with mixed-state statistical models. However, optical flow-based methods, which lend themselves to the extraction of smooth motion fields, and thus might not represent dynamic textures well, which are usually made by chaotic motions in several directions. The time-evolving appearance of textures was explicitly modeled in [START_REF] Soatto | Dynamic textures[END_REF], relying on a Linear Dynamical System (LDS). The method uses the model parameters as the input features. Another well-known family of hand-crafted features are spatio-temporal extensions of the Local Binary Pattern (LBP) method, such as CVLBP [START_REF] Tiwari | Dynamic texture recognition based on completed volume local binary pattern[END_REF] or 2 dimensional LBP computed with 3 orthogonal planes (LBP-TOP) [START_REF] Zhao | Dynamic texture recognition using local binary patterns with an application to facial expressions[END_REF]. [START_REF] Nguyen | Volumes of Blurred-Invariant Gaussians for Dynamic Texture Classification[END_REF] come up with a combination of Gaussian filters and LBP patterns, where LBP descriptors are calculated from both blurred volumes and 3D difference of Gaussians volumes.

Recently, a feature called momental directional pattern (MDP) [START_REF] Nguyen | Momental directional patterns for dynamic texture recognition[END_REF] has been developed.

Since the breakthrough of AlexNet [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], many deep learning approaches have been developed for dynamic texture classification. Some works are based on information that is purely spatial [START_REF] Qi | Dynamic texture and scene classification by transferring deep image features[END_REF], where 2D convolution filters are applied to each frame of a video. In the work of [START_REF] Qi | Dynamic texture and scene classification by transferring deep image features[END_REF], two-level strategy is proposed: utilizing the transfer learning to extract midlevel features and forming a video feature representation by concatenating the mean vector and the diagonal entries of the co-variance matrix of the mid-level features. But such method neglects temporal regularity. In [START_REF] Andrearczyk | Convolutional neural network on three orthogonal planes for dynamic texture classification[END_REF], feature extraction on 3 orthogonal planes based on convolutional neural networks is used. It achieves good results on many dynamic texture benchmark datasets. [START_REF] Hadji | A spatiotemporal oriented energy network for dynamic texture recognition[END_REF] introduce a learning-free ConvNet, operating on oriented 3 rd order Gaussian-based filtering, to extract features, fed to a classifier.

Beside machine learning approaches, data mining techniques, especially sequential pattern mining methods, have also been used to tackle various computer vision problems, such as image classification [START_REF] Nowozin | Weighted substructure mining for image analysis[END_REF]- [START_REF] Yuan | Mining discriminative co-occurrence patterns for visual recognition[END_REF] or action recognition [START_REF] Gilbert | Fast realistic multi-action recognition using mined dense spatio-temporal features[END_REF]- [START_REF] Wang | Mining 3d key-pose-motifs for action recognition[END_REF]. For example in [START_REF] Nowozin | Weighted substructure mining for image analysis[END_REF], a data mining algorithm extracts frequent item patterns, which are used to create a bag-ofvisual-words representation [START_REF] Fernando | Mining mid-level features for image classification[END_REF], for both unsupervised image ranking and supervised image classification. In [START_REF] Fernando | Mining mid-level features for image classification[END_REF], frequent patterns are used as mid-level features for image classification. For video data, [START_REF] Nowozin | Discriminative subsequence mining for action classification[END_REF] propose a pattern-growth mining method to extract sequential patterns. These patterns are then used as a bag-of-words to construct feature histograms, to classify human actions. The presented methods use a hard assignment from input observations to items and patterns. This can lead to unstable patterns and sensitive decisions, which is not suited to the fuzzy nature of dynamic textures. In order to avoid such problems, [START_REF] Wang | Mining 3d key-pose-motifs for action recognition[END_REF] propose a framework for action recognition where patterns are sequences of soft-assignments, rather than items themselves.

In this paper, we propose a novel representation of dynamic texture data, with video classification as final purpose, based on data mining. The proposed framework uses a sequential pattern mining algorithm to discover frequent sequences of image patches, that are referred to as key motifs, inspired from [START_REF] Wang | Mining 3d key-pose-motifs for action recognition[END_REF]. In the training phase, first, a set of symbols, corresponding to representative image patches, is extracted by unsupervised clustering. These symbols are the main modes of patch distribution. Then, for each video class, a set of key motifs is constructed. This extraction handles both deterministic sequences (hard-assigned) and probabilistic sequences (soft-assigned). Sets of key motifs of the different classes are merged into a global one, in which each class is fairly represented. In the inference phase, probabilistic soft-assigned sequences are extracted from the video to be classified. The match between these probabilistic sequences and the key motifs of the unified set, makes a feature vector fed to the final classifier. The experimental analysis studies the use of this motif-based representation for video classification with an extensive parameter study, considering the impact on the number of key-motifs and the impact on classification performance.

We emphasize the fact that the proposed method does not follow the trend of convolutional neural networks, and more generally deep learning approaches. Representing dynamic textures using deep learning can be done by considering the learned features that travel along the successive layers of a convolutional and/or recurrent neural network. Since deep nets achieve very good performances, these features maps do represent well the spatio-temporal patterns that occur in dynamic textures. However, not every deep learning feature is discriminant. Most of deep nets have important redundancy, and interpretating features is intermediate layers is excessively hard, at least because of the quantity. We believe that alternative ways of extracting and representing discriminant image/video patterns, by means of sequential data mining here, should be explored in order to head towards explainable and interpretable machine learning. The code is publicly available at: https://github.com/lphatnguyen/proba seq mining.

II. PROPOSED METHOD

Consider the task of video classification in C classes. In this section, we introduce an efficient framework which is based on a novel mining method of p-sequences, which are sequences of probability vectors. The proposed method is divided into three main stages:

• Unsupervised clustering of patches using a Gaussian mixture model, and construction of p-sequences. • Mining key motifs from p-sequences for each class, thanks to sequential pattern mining. • Constructing feature vectors which are based on probabilistic supports of the union of key motifs and classification using a SVM with χ 2 kernel.

A. Patch extraction and Gaussian mixture clustering

In this part of the paper, a simple method for patch extraction is proposed. From a video V ∈ R T ×H×W , multiple non-overlapping patches of size σ × σ are extracted. From the training dataset, we build a set X = {x i } i=1,...,|X | , where

x i ∈ R σ 2
is a flat vector representation of the i th patch. Theses patches are used for calculating clusters. Beside K-Means clustering, Gaussian mixture models (GMM) are often used as clustering methods in computer vision to create dictionaries of visual symbols. Parameters of the GMM are estimated using the Expectation -Maximization (EM) method. Once they are set, the probability of a given patch can be easily calculated.

By definition, a Gaussian mixture (GM) is a combination of a finite number of K Gaussian distributions N (x|µ k , Σ k ) where k = 1, 2, ...K [START_REF] Nielsen | K-mle: A fast algorithm for learning statistical mixture models[END_REF]. The Gaussian mixture model is made up by mixture weights

π k ∈ R, means µ k ∈ R σ 2 and covariances Σ k ∈ R σ 2 ×σ 2
. The Probability Density Function (PDF) at point x is:

p(x|θ) = Σ K k=1 π k N (x|µ k , Σ k ), 0 ≤ π k ≤ 1, Σ K k=1 π k = 1
where θ is the collection of all parameters of the GM (mixture weights, means and covariances). The expectation maximization algorithm (EM) is used as a learning method to update the parameters in the GM model [START_REF] Nielsen | K-mle: A fast algorithm for learning statistical mixture models[END_REF]. The means computed by the EM algorithm form the dictionary of K symbols. The clustering is done on the whole training dataset X , regardless of video classes.

After having learned the parameters of the model, a soft assignment is used rather than a hard assignment. In other words, a data point is represented by a vector of posterior probabilities rather than an assignment to a cluster. For a given data point x, the posterior probability p(k|x), i.e. the probability that x belongs to the k th cluster, is

p k (x) = p(k|x) = π k N (x|µ k , Σ k ) Σ K j=1 π j N (x|µ j , Σ j ) (1)
Using Eq (1), the soft assignment of a patch of size σ × σ is represented by a probability vector p = [p 1 , p 2 , ..., p K ] ∈ [0, 1] K . At a given spatial position, a sequence of consecutive patches forms a spatio-temporal patch of size T × σ × σ. Applying Eq (1) for each 2D patch of that sequence and each cluster k, the p-sequence P of length T is built, such as: P = p 1 , p 2 , ..., p T where each p j (j = 1, 2, ..., T ) is a probability vector in [0, 1] K . Then, P is the set of p-sequences to be mined, P = P 1 , ..., P |P| .

B. Mining key motifs

Motif mining is a method to discover relevant subsequences as patterns in a sequential database. Many early motif mining methods are based on the Apriori algorithm for databases where items are certain. Such algorithm can be considered as a method of pattern expansion. The main idea is that subsequences of frequent motifs can also be frequent. Algorithm 1 demonstrates how super-sequence patterns can be generated from shorter mined patterns. In other words, from patterns containing l symbols, referred to as l-motifs, candidate patterns of l + 1 symbols are generated. For a particular case, an item can be considered as a 1-motif composed of one symbol in the dictionary.

Unlike in itemized sequence pattern mining, our data to be mined (p-sequences) are uncertain (probabilistic). To do so, we follow the mining procedure proposed in [START_REF] Wang | Mining 3d key-pose-motifs for action recognition[END_REF]. Algorithm 2 describes the motif-mining algorithm applied in our case for each class c, where each element t i , i = 1, ..., |T | of the candidate motif T is a symbol corresponding to a cluster k = 1, ..., K of the learned GM model. Notice that the Algorithm 1: Expansion Algorithm

Input: T l = {T l 1 , T l 2 , ..., T l |T l | }: set of l-motifs Output: T l+1 : set of l + 1-motifs 1 T l+1 ←-∅ 2 for i = 1, 2, ..., T l do 3 tail ←-T l i (2 : T l i ) 4 for j = 1, 2, ..., T l do 5 head ←-T l j (1 : T l j -1) 6 if tail = head then 7 T l+1 ←-T l+1 + concat(T l i , T l j ( T l
η(T , P ) = max M |T | i p M (i) t i , Subject to:M (i) < M (i + 1), M (i + 1) -M (i) ≤ g. (2)
where g is the maximum gap constraint and p

M (i) t i
is the probability, at temporal position M (i), of symbol t i . Considering that the candidate motif T has |T | elements and the p-sequence P has T probability vectors, as described in [START_REF] Wang | Mining 3d key-pose-motifs for action recognition[END_REF], the probabilistic support function f (T (1 : |T |) , P (1 : T )) measures the probability of matching between T (1 : |T |) = t i i=1,2,...,|T | and P (1 : T ) = p j j=1,2,...,T . At a given position i in T and j in P , the probabilistic support function f (T (1 : i) , P (1 : j)) is calculated using equation 3.

f (T (1 : i) , P (1 : j)) = p j t (i) × max k∈{j-g,...,j-1} f (T (1 : i -1) , P (1 : k)) (3)
which can be implemented easily using dynamic programming. When finished constructing the dynamic table, the probabilistic support is the maximum value of the last row of the table: η(T , P ) ←-max j∈{1,...,|P |} f (T (1 : |T |) , P (1 : j)) .

C. Classification

A video sequence contains multiple non-overlapping video subsequences of size T × H × W . Each video subsequence consists of several non-overlapping spatio-temporal patches (represented by the P using Eq. 1). In this part, a bag-ofmotifs method is used from the mined motifs of all classes. This process is divided into 3 steps: (a) uniting mined motifs of each class, (b) constructing probability histograms for each video sample and (c) classification learning. After having mined sequential motifs using Algorithm 2, a union of the motifs of all considered classes is made. The union set of mined motifs of all classes is T ∪ = C c=1 T c . The previous mining of key motifs on a per-class basis ensures that, in this union, each class is sufficiently represented. If key motifs had been extracted regardless of classes, we could end up with a set where some classes would be much more represented than other ones when dealing with unbalanced datasets. Moreover, making this union allows us to define a global bag of motifs that will take into account the relative representation of each motif in each class, including the fact that some motifs could be ambiguous (i.e. represent different classes).

A bag-of-feature vector is generally a histogram that counts the number of occurrences of each symbol in a sample [START_REF] Leung | Representing and recognizing the visual appearance of materials using three-dimensional textons[END_REF]- [START_REF] Lo | Classification of software behaviors for failure detection: A discriminative pattern mining approach[END_REF]. However, in our case, a feature vector is built based on the probabilistic supports. For a given test p-sequence P , a feature vector f , is made up of the probabilistic supports between each key motifs of the union set T ∪ and P . Hence, we have f ∈ R |T∪| . To train the SVM classifier and at test time, we use feature vectors of video subsequences, instead of feature vectors of individual p-sequences. A feature vector of a video subsequence is calculated using the average of the f 's of the p-sequences in this video subsequence. For the classification step, an SVM classifier based on the χ 2 kernel is used [START_REF] Wang | Dense trajectories and motion boundary descriptors for action recognition[END_REF].

D. Differences with Wang et al. [19]

There exists several differences between our method and [START_REF] Wang | Mining 3d key-pose-motifs for action recognition[END_REF]. Since the application is different, we use image instead of 3D body joint coordinates. Second, we extract symbols by fitting a GMM, while [START_REF] Wang | Mining 3d key-pose-motifs for action recognition[END_REF] use the activated simplices method [START_REF] Wang | Recognizing actions in 3d using action-snippets and activated simplices[END_REF]. Third, contrary to what is done in [START_REF] Wang | Mining 3d key-pose-motifs for action recognition[END_REF], the mining is done with a fixed threshold ϵ. Indeed, thanks to the modification mentioned above in Algorithm 2, there is no need to change ϵ during the mining process in order to obtain a given number of key-motifs. Experiments show that this threshold does not play a crucial role for accuracy but only for compactness of the model (number of extracted motifs). Lastly, [START_REF] Wang | Mining 3d key-pose-motifs for action recognition[END_REF] classify a test example by assigning it to the class which has the key motifs leading to the highest probabilistic supports. We found that characterizing a test example with the probabilistic support of all key motifs to be a more discriminant representation for classification.

E. Additional study on FFT magnitude patches

In our research, we are interested in exploring the magnitude values of the Fourier transforms of the texture input patch instead of the gray-scale patches. In other words, the Fourier transform is applied for a single patch but not the entire frame. The final Fourier transform magnitude patch x F T is computed using the following procedure: (1) each gray-scale patch is transformed into the frequency domain using the FFT algorithm which is widely used for computing the Fourier transform of a patch x F T ∈ C σ×σ ; (2) a top left quadrant of the FFT patch is taken as it is symmetric along the horizontal and vertical axis, so the FFT patch has a size of C σ 2 × σ 2 ; (3) the magnitude of the FFT patch is computed as

x F T ∈ R σ 2 × σ 2 .
Then, the key motif mining approach is used similarly on FFT patches instead of gray-scale patches.

III. EXPERIMENTS

In this section, we are presenting the results obtained by our proposed method on different benchmark datasets (Traffic, UCLA) while comparing them to those of state-of-the-art methods.

A. Datasets and protocols

Traffic dataset [START_REF] Chan | Probabilistic kernels for the classification of auto-regressive visual processes[END_REF] consists of 254 video samples. The sequences are recorded with a resolution of 320 × 240 with a temporal length between 42 and 52 frames, at 10 fps. Each video sample is then spatially resized to 80 × 60 and cropped to a size of 48×48 over the area where the motion is the most prominent. There are 3 classes in this database: heavy, medium and light traffic. For this dataset, a 4-fold evaluation protocol is used [START_REF] Chan | Probabilistic kernels for the classification of auto-regressive visual processes[END_REF]. The average score of the 4-fold is recorded as the final result.

UCLA dataset [START_REF] Soatto | Dynamic textures[END_REF] consists of 50 classes, where each class contains 4 sequences. As the result, the dataset has 200 dynamic textures sequences in total. Each original sample is captured with 75 frames, each of size 160 × 110. A slightly modified version of UCLA dataset is often utilized for dynamic texture classification where the original samples into sub-sequences of size 48 × 48. Three popular challenges of this dataset [START_REF] Nguyen | Momental directional patterns for dynamic texture recognition[END_REF], [START_REF] Tiwari | Dynamic texture recognition using multiresolution edge-weighted local structure pattern[END_REF], [START_REF] Ribas | Dynamic texture analysis with diffusion in networks[END_REF] are often used for classification. In the 50-class (4-fold) configuration, a quarter of the data in each class is addressed as testing set and the remaining for the learning. The experiment is repeated 4 times and the average score of the 4 folds is reported as the final result. As for 9-class configuration, a reorganization of the original 50-class configuration is made into 9 classes of boiling water, fire, flower, fountains, plants, sea, smoke, water and waterfall. Half the videos in each class are randomly selected for training and the remaining for testing. This step is repeated 20 times and the mean score of all the repetitions is recorded as the final result. In the 8-class configuration, following the 9-class configuration, the plants class is omitted since the number of samples is much greater than those of other classes. The same evaluation metric as the 9-class configuration is used in order to evaluate this sub-dataset.

B. Experimental setup

As mentioned above, depending on the dataset and the configuration, a N -fold validation is performed. Each fold is composed of different set of videos. In our experiment, a video sample of each dataset is divided into multiple nonoverlapping sequences of temporal lengths of 1 second. For the traffic database, a video subsequence contains 10 frames. For the UCLA dataset, a video subsequence has 15 frames. The final classification result for each test video sample is the majority vote of classified subsequences in the video sequence, which means the final result of a video is the majority class assigned to the subsequences of the video sample.

For the learning of the GM model, 4 patches in each video frame from training set are randomly selected. The GM model is initialized randomly. To assess the stability of our method against the choice of symbols (that depends on these random selections), the impact on the results is evaluated over 8 runs for each split of the N -fold validation.

The minimum support threshold ϵ (algorithm 2) is an important parameter. Therefore, an analysis could be useful to study the impact of such parameter ϵ on the number of mined motifs. Consequently, during the 8 runs, an analysis of parameters is also conducted to evaluate the impact on both accuracy and the number of motifs. This process is done for each set of parameters on the Traffic dataset and UCLA dataset. This analysis is only done on the gray-scale patches.

In the training and testing stage, the spatio-temporal patches are extracted at every possible position with a non overlapping condition. For each dataset and their configurations, both grayscale patches and FFT patches are evaluated.

C. Experimental results

Traffic database: Looking at figure 1, the best results are obtained with a patch size of 8 (compared to 12). This is the most important parameter. Figure 1 shows that the number of motifs generated from codebook of sizes 20 and 50 are almost the same. Furthermore, the performances based on these codebook sizes are very close. For the minimum support threshold, good results can be achieved with a minimum support between 0.01 and 0.07. However, the number of mined key-motifs increases as the minimum support threshold decreases. Therefore, rather than using 500 or 600 key motifs with ϵ of 0.01, high results (eg. 96,56% of accuracy) can be achieved by using approximately 100 key motifs with an ϵ = 0.05. Column 5 of table II shows the performance of Datasets UCLA Traffic Method 50-4fold 9-class 8-class KDT-MD [START_REF] Chan | Classifying video with kernel dynamic textures[END_REF] 97.5 ---DFS [START_REF] Xu | Dynamic texture classification using dynamic fractal analysis[END_REF] 89.5 ---CFV [START_REF] Wang | Chaotic features for dynamic textures recognition[END_REF] -85.10 85.00 -3D-OFT [START_REF] Xu | Scale-space texture description on sift-like textons[END_REF] 87.10 97.23 99.50 -HLBP [START_REF] Tiwari | A novel scheme based on local binary pattern for dynamic texture recognition[END_REF] 95.00 98.30 97.50 -MEWLSP [START_REF] Tiwari | Dynamic texture recognition using multiresolution edge-weighted local structure pattern[END_REF] ' * ' indicates that the result is averaged over 8 runs.

the proposed method for the Traffic database. The proposed method outperforms most of existing methods on this dataset and it is even on par with the diffusion-based method [START_REF] Ribas | Dynamic texture analysis with diffusion in networks[END_REF] with only 0.04% less (96.60% and 96.56%) with the average accuracy over 8 runs. Comparing to the diffusion-based model methods [START_REF] Gonc | Dynamic texture analysis and segmentation using deterministic partially self-Avoiding walks[END_REF], [START_REF] Gonc ¸alves | A complex network approach for dynamic texture recognition[END_REF], our motif-based method shows better results by at least 1%. Looking at the results obtained by the LBP-based methods [START_REF] Zhao | Dynamic texture recognition using local binary patterns with an application to facial expressions[END_REF], [START_REF] Zhao | Dynamic texture recognition using volume local binary patterns[END_REF], it classifies the dataset much better by up to almost 3% being more compact in terms of feature vector size, which can be seen with Table I (16384 [START_REF] Zhao | Dynamic texture recognition using volume local binary patterns[END_REF] and 768 [START_REF] Zhao | Dynamic texture recognition using local binary patterns with an application to facial expressions[END_REF] to 118 with our approach). In addition, the result obtained by using a quadrant of FFT image as input is very high with 98.44% of accuracy (table II). It is currently the highest score for this dataset. This result shows the efficacy of the proposed method for the Traffic dataset.

UCLA 50-class: Best results are obtained with the patch size of 8 (figure 2). But this time, the accuracy is better with 50 clusters than with 20 clusters. This can be explained by the complexity of textures occurring in the videos. For the minimum support, good results are achieved with ϵ between 0.19 and 0.3. In figure 2, the averaged highest result is obtained with only 651 motifs. From table II, the proposed method achieves 98.50% of the classification rate on this subset as the average score over 8 runs. It performs on par with the diffusion-based model (DM) [START_REF] Ribas | Dynamic texture analysis with diffusion in networks[END_REF] and outperforms other recent existing methods for example HLBP [START_REF] Tiwari | Dynamic texture recognition based on completed volume local binary pattern[END_REF] or MEWLSP [START_REF] Tiwari | Dynamic texture recognition using multiresolution edge-weighted local structure pattern[END_REF] as well as many other hand-craft methods. Nonetheless, in comparison with the DL method, the mining of key motifs method is just 1% less than the existing DT-CNN models [START_REF] Andrearczyk | Convolutional neural network on three orthogonal planes for dynamic texture classification[END_REF]. Furthermore, in comparison to a recent LBP-based method, MDP [START_REF] Nguyen | Momental directional patterns for dynamic texture recognition[END_REF] which scores 100% of accuracy for this subset, our method is just 1.5% less, with a higher compacity from 3880 to 651 features ( see Table I). When FFT patches are used, the result (99.50%) increases by 1% in terms of accuracy and even on par with the DL approaches.

UCLA 9-class and 8-class: Figures 3 and4 illustrate the parameters analysis of the UCLA 9-class and UCLA 8-class schemes. The tendency of classification rates is approximately the same as the UCLA 50-class and Traffic dataset as good results are obtained with ϵ ranging in the middle of the analyzing interval. The number of motifs begins to converge as ϵ passes 0.15. From table II, our model outperforms some of the other classical computer vision methods of about 1% or even better than CFV [START_REF] Wang | Chaotic features for dynamic textures recognition[END_REF] of about 10% for the 9-class configuration. Nevertheless, it cannot achieve as good results as DL models [START_REF] Andrearczyk | Convolutional neural network on three orthogonal planes for dynamic texture classification[END_REF] or a recent vision-based model of MDP [START_REF] Nguyen | Momental directional patterns for dynamic texture recognition[END_REF] which are at 98.35% and 98.90% respectively. In the 8class configuration, our method reaches a result of 95.45% of accuracy. This result is higher than many methods such as CDT-TOP [START_REF] Gonc | Dynamic texture analysis and segmentation using deterministic partially self-Avoiding walks[END_REF], CVLBP [START_REF] Tiwari | Dynamic texture recognition based on completed volume local binary pattern[END_REF] or CFV [START_REF] Wang | Chaotic features for dynamic textures recognition[END_REF]. Nonetheless, the best results for this configuration is obtained with GoogleNet [START_REF] Andrearczyk | Convolutional neural network on three orthogonal planes for dynamic texture classification[END_REF] with an accuracy of almost 100%. Plus, using only one quadrant for this scheme also results in good result (96.55%) but not as high as with grayscale images as inputs. Moreover, UCLA 8-class with FFT patches helps improving the performance by only 0.1% (95.54% in accuracy) comparing with the grayscale patches. Overall, the results, even if not the best, are satisfactory considering both accuracy and compacity (Table I).

In general, the performance of the proposed approach is comparable to the state-of-the-arts methods. Looking at the overall accuracy of the experiments, different sets of parameters are tested in order to measure the impact of each parameter for each module. Overall, it can be seen that the two most important parameters in the approach is the patch size σ and the number of cluster K. A smaller patch size of 8 with a greater number of cluster of 50 seems to perform best. Moreover, a maximum gap constraint also contributed to the performance of the approach. A smaller gap seems to generate less motifs than a bigger gap (a maximum gap constraint of 3 and 5) but with slightly better accuracy. This means that the set of extracted key motifs using a gap of 3 are more meaningful and more discriminant than with a gap of 5. As for the minimum support threshold (directly related to the number of motifs), increasing gradually this value can help us reach an optimal value and achieve highest result. This helps to remove some redundant motifs extracted from the mining process and to keep only the meaningful motifs. Passing this optimal value, the system begins to remove some discriminant motifs between classes and the accuracy decreases.

IV. DISCUSSION AND FUTURE WORKS

Experimental classification results show that our key motifsbased representation is relevant for dynamic textures. The proposed method performs on par existing work, including deep learning approaches.

For further perspectives, mining motifs using p-sequences can be effective but not all mined motifs are useful for visual representation of video data. Therefore, motif filtering could be applied in order to extract only useful motifs. This addition should help speeding up the classification step and reduce the number of elements in the feature vector. One way to do it is to compute an entropy-like value of each motif. This value can be computed using the probabilistic support computed using p-sequences in each class and the whole dataset. If the entropy of the motif is low, the extracted motif is meaningful and represents well the input video. On the other hand, if the computed motif entropy is high, the motif appears in every class with an equal distribution to the data in every class. Such motifs with high entropy should be eliminated from the motif set.

Algorithm 2 : 4 T 6 support ←-0 7 for j = 1 , 2 ,

 246712 Key motifs Mining Algorithm Input: P c : p-sequences set of class c, ϵ: support threshold Output: Mined key motifs of the class c-T c 1 T 1 ←-{1 -motifs} 2 l ←-2 3 while T l-1 ̸ = ∅ do l ←expand(T l-1 ) 5 for i = 1, 2, ..., T l do ..., |P c | do 8 support ←support + η T l i ,

12 Fig. 1 .

 121 Fig.1. An analysis of the stability of the method in terms of minimum support threshold and the number of mined motifs for the traffic dataset.

  motifs gap = 3, K = 20, σ = 8 gap = 3, K = 50, σ = 8 gap = 3, K = 20, σ = 12 gap = 3, K = 50, σ = 12 gap = 5, K = 20, σ = 8 gap = 5, K = 50, σ = 8 gap = 5, K = 20, σ = 12 gap = 5, K = 50, σ = 12

Fig. 2 .

 2 Fig.2. An analysis of the stability of the method in terms of minimum support threshold and the number of mined motifs for the UCLA 50-class dataset.

  motifs gap = 3, K = 20, σ = 8 gap = 3, K = 50, σ = 8 gap = 3, K = 20, σ = 12 gap = 3, K = 50, σ = 12 gap = 5, K = 20, σ = 8 gap = 5, K = 50, σ = 8 gap = 5, K = 20, σ = 12 gap = 5, K = 50, σ = 12

Fig. 3 .

 3 Fig.3. An analysis of the stability of the method in terms of minimum support threshold and the number of mined motifs for the UCLA 9-class dataset.

Fig. 4 .

 4 Fig.[START_REF] Lu | Dynamic texture recognition by spatio-temporal multiresolution histograms[END_REF]. An analysis of the stability of the method in terms of minimum support threshold and the number of mined motifs for the UCLA 8-class dataset.

  Algorithm 1) is called at line 4. Algorithm 2 depends on the computation of the probabilistic support η (T , P ). This support measures how well the motif T is matched to a p-sequence P . The calculation of the probabilistic support searches for a mapping table M where each M (i) ∈ {1, 2, ..., |P |} represents a row i of the computed dynamic table. The probabilistic support is computed as

		j )
	8	end
	9	end
	10 end
	expansion step (

TABLE II COMPARISON

 II OF RECOGNITION RATES (%) ON THE UCLA DATASET AND TRAFFIC DATABASE.


Jose Alfredo Guerrero
email: jguerrer@ieee.org

Juan Antonio Escareño

Yasmina Bestaoui

I A Escareno

Quad-rotor MAV trajectory planning in wind fields

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

Quad-rotor MAV Trajectory Planning in Wind Fields

1.A. Guerrero, 1.A. Escareno, Y. Bestaoui

Ahstract-This paper addresses the problem of time optimal path planning for a quadrotor helicopter evolving in a region of known winds. Usually, the flight control of quad rotors subject to wind disturbances challenge seeks to find the optimal control to keep track of a desired trajectory in a windy region. This approach has one major disadvantage: the quadrotor flight control has to compensate for trajectory deviations; therefore, the energy consumption becomes an issue. Most unmanned aerial vehicles (UAV) navigation techniques use waypoints to accomplish their missions. In the framework of a waypoint based navigation, a promising path planning strategy would be a time optimal approach in which the UAV would take advantage of wind to reach its next waypoint; therefore, saving time and energy (under constant forward velocity constraint).

A model separation is used to simplify the control of the six degrees-of-freedom (6DOF) dynamics of the quadrotor. Such approach allows to deal with quadrotor's 3D-motion through two subsystems: dynamic (altitude and MAV-relative forward velocity) and kinematic (nonholonomic navigation) subsystems.

In terms of control, a hierarchical control scheme is used to stabilize dynamic and kinematic underactuated subsystems in volved in the navigation task. The time optimal path planning is computed using a dynamic optimization method for continuous systems with some state variables specified at an unspecified terminal time. Results have been validated in simulation.

I. INT RO DUCTION

The applications of Miniature Air Vehicles (MAVs) have widely diversified during the last years. They comprise both military and civilian, though the latter has had a lower development rate. The application scenarios range from homeland security, pipeline surveillance, nuclear facilities monitoring, natural disaster damage assessment, etc. These real-world missions evolve in outdoors environments where the autonomous aerial vehicles are exposed to adverse atmo spheric conditions; most of times, a windy environment. The majority of real-world MAV missions evolve in a waypoint navigation framework. An intuitive solution is to use a wind compensation approach which aims to compensate MAV trajectory for disturbances due to meteorological phenomena; therefore, energy consumption becomes an issue. For this reason, navigation methods that take advantage of wind velocity and direction are required to achieve the guidance objective between two any given waypoints; therefore, saving energy.

The literature in MAV guidance in regions of mild winds is vast and particularly focused on fixed-wing configurations and blimps. As in [1] where the authors use a linear quadratic regulator for optimal guidance of a fixed-wing MAY. The controller, evaluated at simulation level, aims to follow a straight-line and circular paths under different constant wind values. In [2] is experimentally implemented a path follow ing controller using constant-velocity MAV under moderate wind conditions (20% to 50% of MAV airspeed). In [3] is presented a path planning and control algorithms meant to survey multiple-waypoints while considering heading and constrained rate heading. The path planning is treated as an optimization problem and assumes the knowledge of the constant wind component, while the sliding surface controller is aimed at dealing with small time varying wind compo nents. A nonlinear controller is presented in [START_REF] Brezoescu | Straight-line Path Following in Windy Conditions[END_REF] regarding to the following of a straight-line reference considering various constant lateral wind components (crosswind) up to 40% of the MAV speed.

Time optimal aircraft guidance in presence of winds has been discussed in [START_REF] Bijlsma | Optimal Aircraft Routing in General Wind Fields[END_REF] and [START_REF] Jardin | Neighboring Optimal Aircraft Guidance in Winds[END_REF]. [START_REF] Bijlsma | Optimal Aircraft Routing in General Wind Fields[END_REF] presents an aircraft routing algorithm to determine a minimal time path for an airplane using a Bolza problem approach. In [START_REF] Jardin | Neighboring Optimal Aircraft Guidance in Winds[END_REF], a neighboring optimal control approach has been adopted to obtain near optimal trajectories using nominal solutions to the Zermelo problem. Both approaches consider the kinematic model of the aircraft, evolving in the 2D space, in the optimization problem. Optimal path planning between two waypoints considering the presence of obstacles and no-fly zones have been discussed in [START_REF] Mattei | Smooth Flight Trajectory Planning in the Presence of No-Fly Zones and Obstacles[END_REF] and [START_REF] Prats | Lexicographic Opti mization for Optimal Departure Aircraft Trajectories[END_REF]. In [START_REF] Mattei | Smooth Flight Trajectory Planning in the Presence of No-Fly Zones and Obstacles[END_REF] a trajectory planning is presented. This approach considers the use of core path graphs and optimal control in order to obtain an optimal tra jectory between the origin and the destination while avoiding both obstacles and no-fly zones. The results were validated in simulation considering 2D and 3D cases. A lexicographic approach to optimal trajectory planning for the departure of aircrafts considering no-fly zones has been addressed in [START_REF] Prats | Lexicographic Opti mization for Optimal Departure Aircraft Trajectories[END_REF]. By using this approach, the optimization objectives can be ordered hierarchically. Recently, in [10] a trajectory planning for bridge inspection using harmonic potential functions has been presented. [11] presents an analysis of the time optimal trajectories of an airship. The time optimal problem is formulated as a generalization of the Zermelo's navigation problem. An analysis of the Lagrange multipliers and the set of solutions is presented.

This paper presents a hierarchical control architecture for time optimal quadrotor helicopter flying in a region of a-priori known winds. A simplified 3D model is used by splitting the 6DOF dynamics into two subsystems: a dynamic (forward velocity and altitude) subsystem and a kinematic (navigation) subsystem. Likewise, a hierarchical control scheme is used to stabilize underactuated dynamics 1 and kinematics involved in the navigation control task. A nonlinear bounded controller is used to stabilize the angular error dynamics regarding to avoid singularities and actuators limits. A time optimal path planning is developed considering an evolution of the aerial robot in the 3D space. The resulting path is then used as the reference path to be followed by the nonlinear bounded controller. Simulations results shows the validity of the proposed approach.

This work is organized as follows: the dynamic model of the quadrotor is presented in section 2. The time-optimal path planning is presented in section 3. Section 4 presents the control strategy. Numerical simulations results are presented in Section 5. Conclusions and perspectives are finally given in Section 6.

II. QUA DROTO R DYNAMIC MO DEL

A. Quadrotor Kynodynamic Model

In this section a two-rotational coupled model is presented to drive the aerial robot through the 3D space, with the aim of simplifying the navigation control task performed by the vehicle (Fig. 1).

• A regulated roll motion having small angular fluctua tions is assumed.

• Pitch rotation n O defines a dynamic subsystem contain ing altitude and MAV 's relative motions.

• Yaw rotation n 'lj; applied to velocity vector V r yields a kinematic motion subsystem. where Vr stands for the magnitude of the linear velocity relative to the quadrotor V n T is the magnitude of Tb, m is the mass of the vehicle which in the following is normalized to 1, T O is the pitch control input.

2) Kinematic subsystem (K -plane): The transformation n 'lj; , about et , drives the dynamic plane according to head ing motion. Therefore, T� generates a planar translational motion along the inertial axes exe y , called from now on K plane (Kinematic plane). Provided that gravity is treated in the dynamic plane, then, it is more convenient to address vehicle's motion from a kinematic viewpoint, i. e. focusing on motion subsystem resulting from heading the relative velocity vector V r. The latter is expressed by the following (2) T 'lj; where T 'lj; is the yaw control input.

In the following section, the time optimal path planning for an aerial vehicle evolving in the 3D space is presented. It will be shown that the kinematic model of the quadrotor helicopter can be obtained from the general 3D kinematic model of an aerial vehicle.

III. PATH PLANNING IN WIND FIELDS

Lets consider the problem of a quadrotor helicopter flying in region of known winds available through meteorological forecasts or doppler radar measurements. The magnitude and the direction of the winds are known to be functions of position, i.e. Wx = Wx(x, y, z), Wy = Wy(x, y, z) and Wz = Wz(x,y,z) where (x,y,z) are rectangular coordi nates and (Wx, Wy, Wz) are the velocity components of the wind. The quadrotor velocity relative to the air mass is V, a constant. The goal is to find the minimum-time path between a pair of waypoints. The kinematic model of an aerial vehicle is given by :i;

if i V cosxcos, + Wx, V sin X cos, + Wy, Vsin,+ Wz (3)
where X is the heading angle of the aerial vehicle relative to the inertial frame, , is the flight path angle. Notice that the kinematic model of the quadrotor in (2) previous section can be obtained from (3) using, = O. The heading angle X is approximated by the yaw angle 1jJ and the velocity relative to the air mass V is approximated by the quadrotor forward velocity -Vr.

The time optimal trajectory generation can be formulated as follows: subject to

:i; if i Ul(t)+WX, U2(t) + Wy, U3(t) + Wz (4) (5)
where Ul (t) = V cos X cos" U2(t) = V sin X cos, and U3 (t) = V sin,; with the constraint ui(t) + u §( t) + u�(t) ::; V';a x . If the final waypoint is reachable at any time, then it is reachable in minimal time. However, if W » V, i.e. the wind is too strong with respect the quadrotor velocity, there may be points that are not reachable at all. The Hamiltonian of the system is H =Ax( V cos X cos /' + Wx) + Ay( V sin X cos /' + Wy) +Az( V sin /' + Wz) + 1. [START_REF] Bijlsma | Optimal Aircraft Routing in General Wind Fields[END_REF] where A x, A y and A z are the Lagrange multipliers.

The Euler-Lagrange equations are ,\ __ oH _ -A oWx _ A oWy _ A oWz

(lO) (11)

After some mathematical manipulation, it is easy to show

A. Constant Wind

From (15) note that, if the Wx, Wy and Wz are constant, implies that X = const., i.e., the minimum-time path are straight lines.

B. Linear Variation of Wind Velocity

Let Vw be the wind velocity, a constant. If Wx(x, y) =t=VwY and Wy(x, y) = 0, it has been proved in [START_REF] Jardin | Neighboring Optimal Aircraft Guidance in Winds[END_REF] that

y = � Ci �X -sin\j ) , (17) x = � ( asinh(tanXj ) -asinh(tanx) +tanx ( _ . _ 1 __ _ . _ 1_ ) 2Vw smx smXj _ _ 1 _ (tanXf -tanx) ) . (

18) COSXj

where Vw is a positive constant, the wind velocity magnitude. The evolution of the heading angle X and flight path angle /' is reduced to X = ±Vw cos 2 ( X) ' -cos X cos/'

(12) A A y -sin X cos/' (13) A A _ -sin/, z - A . ( 14
)
where A = V + WxCOSXcos/,+ Wysinxcos/,+ Wz sin/,.

Introducing (12), ( 13), ( 14) into (7), ( 8), ( 9), we obtain ,\ _ cos X cos /' oWx sin X cos /' oWy sin /' oWz 'Y = ± Vw sin /' sin X cos /' cos X.

(

) 20
The time to go is given by

1 T = V, (tanXf -tan X). w (21)
Notice that in the case of a quadrotor helicopter, the heading angle will be used as reference input for the kine matic control subsystem discussed in the following section.

It is worth to mention that other vehicles may use both heading angle and flight path angle as reference inputs for the trajectory tracking control.

IV . CONTROL STRATEGY

The navigation control design relays on three modules or subsystems: D-plane, K-plane and path planning sub systems. Within D-plane the controller aims to define the relative MAV velocity vector V r and desired altitude, while in the K -plane whose controller aims to steer V r so that the vehicle follows the desired path given by the path planning module. These subsystems are coupled via V r generated by the dynamic subsystem and fed to the kinematic subsystem to perform the navigation task and 'ijJd generated by the path planning module and provided to the kinematic subsystem as a reference input (Fig. 2).

A. Dynamic controller

At this level, our aim is to synthesize a control algorithm to

The controller used to attain the path-following objective follows a similar structure of (35), i.e.

T'l/J = -kv,p �k p ,p tanh � + 1/; d

As a result, it is guaranteed that the state vector (�,�)

converges asymptotically to the origin, using similar stability analysis used for e. The it is possible to guarantee that the flying robot reaches the desired path.

V. SIMULATION RESULTS

To illustrate the proposed solution, numerical simulations were performed to evaluate the quad-rotor flight performance while navigating through a region of winds. In the case, the initial velocity of the quadrotor is Vr = 0.4m.!sec, the initial altitude is O.gm., the linear varying wind velocity is 0.022m1s. The desired altitude was fixed to 1m. while the The dynamic controller uses vr as desired the velocity Concerning the kinematics evolving within K-plane, the hierarchical controller regulating the path deviation kine matics via the heading dynamics presents a satisfactory performance. Figure 7 illustrates the kinematic controller. an impact on the tracking error.

VI. CONCLUDING REMARKS AND FUTURE WORKS

The present paper proposes a model reduction regard

1)

  Dynamic subsystem ( D-plane): This motion subsystem arises from the coordinate transformation n o , about e2, applied to the main thrust vector Tb, which leads to the well known 3DOF longitudinal dynamics. Translational and rota tional dynamics evolve within the plane resulting plane e� e� called from now on D-plane. The advantage of this plane is that gravity might be included in the dynamic equations since the coordinates transformation from the inertial to the actual frame has no effect on gravity vector. The corresponding equations are given next -T sin B -TcosB + mg (1) T O

Fig. 1 .

 1 Fig. 1. D-K Planes

  __ oH _ -A oWx _ A oWy _

  that A x, A y and A z are given by and A x

+++

  A ox',\ _ cos X cos /' oWx sin X cos /' oWy sin /' A oy' ,\ _ cos X cos /' oWx sin X cos /' oWy sin /' oWz z -A oz• From previous equations it is easy to obtain an expression that describes the evolution of the heading angle and flight + sin /' sec /' sin X ox -cos X 8y -cos X 8y' 'Y cos X--+ CosXCos'Ysm'Ysmx-ox oy 2 oWy 2 oWz + sin X cos 'Y sin 'Y --+ sin 'Y sin X-oy oy -cosxcos 'Y--+ sinxcos 'Y--+ sin 'Ycos'Y--. oz oz oz (16)

  achieve the desired references of drone-relative velocity and altitude, that is to say, errors Vq = Vr -v� and z = Z -zd are rendered to zero. The desired forward velocity and altitude are defined as Vd ' control of the quad-rotor kinematics, this introduces the following error variable (40) with 1/J d given by the path planning subsystem discussed in section III. Notice that the kinematic model of the heading subsystem is given by 1/; d = =t=Vwcos 2 1/J (41) The control of the heading error dynamics (1[;) requires the time derivative, acceleration, of (41) to guarantee the trajectory following. The corresponding error dynamics is thus written as

  desired cruise velocity was set to 0.5m.!sec. Extensive cases have been studied, however, due to space limitations, three trajectories departing from the same coordinates (-4, 4) are discussed in this section. The final waypoints are (40, 4) , (40, 0) and (40, -4) . Using the Zermelo-based approach the initial and final headings as well as the time to go from departure waypoint to each final waypoint is presented in table I. The corresponding trajectories are depicted in Figure 3.

  reference for the outer-loop controller. The underactuated subsystem vr-B is successfully controlled through a hier archical control scheme. The quad-rotor altitude controller achieves a constant desired altitude of 1m. The corresponding results are shown on Figures 4-6.

Fig 5 xlmlFig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .

 54567 Fig. 3.

Fig. 8 .

 8 Fig. 8. Trajectory tracking in a region of winds with destination (40,4).

  ing to simplifying the control task by partitioning the the complete dynamic model into a dynamic subsystem (D plane) that focus to control forward velocity and altitude, and a kinematic subsystem (K-plane) that allows to treat the quad-rotor navigation problem as a nonholonomic-like path following. From this model simplification arises the possibility of applying a zermelo-based navigation scheme in the quad-rotor navigation flight with satisfactory results in simulations. The following objective is to validate the proposed approach by means of an experimental air plat form to perform autonomous navigation in wind controlled situations. It is important that an outdoors application of the approach is considered by the authors, which implies an improvement in the kinematic subsystem control by designing robust algorithms that accounts the presence of external disturbances (e.g. wind gust).VII. ACKNOW LE DGEMENTSJ.A. Guerrero thanks IBISC laboratory, Universite d'Evrywhere he developed the time-optimal path planning for UAVs. He also thanks the LISTIC laboratory, Poly tech Annecy for the support to complete this work. RE FE RENCES [I] A. Ratnoo, P.B. Sujit and M. Kothari, "Adaptive Optimal Path Fol lowing for High Wind Flights", IFAC world Congress, Milano, 2011. [2] D.R. Nelson, D. Blake, T Barber, W. McLain, R.w. Beard, "Vector Field Path Following for Small Unmanned Air Vehicles", IEEE American Control Conference, Minneapolis, 2006. [3] TG. McGee, J.K. Hedrick, "Path Planning and Control for Multiple Point Surveillance by an Unmanned Aircraft in Wind", IEEE American Control Conference, Minneapolis, 2006

  :------: :------: :------: , 7 0 ---: ' ':-5 --:2':-0 --:2':-5 --:3':-0 --:3 ':-5 ---:' ':-0 -�'•5 '1m]

Fig. 9 .Fig. 10 .

 910 Fig. 9. Trajectory tracking in a region of winds with destination (40,0).

TABLE I

 I TIME TO GO. INITIAL AND FINAL HEADINGS DEPARTING FROM (-4,4)

	Destination (m) (40,4)	Time (sec) 90.5505	XO (deg) 315.13	Xj (deg) 44.87
	(40,0)	84.8267	313.29	38.76
	(40,-4)	80.5021	310.93	31.64

D-plane

where ar and az are positive constants. To stabilize the sub system (1) we assign a linear stable behavior to translational dynamics �DT through polynomials '1(Vr) and 1 2 ( 2,i) as presented in [START_REF] Sepulchre | Constructive Nonlinear Contro[END_REF].

-Tsine -Tcose + g

Such stable behavior is achieved provided that desired values are achieved

where, '1 = k P r vr and '2 = k p z 2 + kvz i, with k P r' k p z To design the control law for the inner-loop heading dynam ics, we proceed similarly as in [START_REF] Ngo | Integrator backstepping design for motion systems with velocity constraint[END_REF], where a backstepping based controller is obtained considering a bounded velocity state (see A2 and A3).

.

now, using sech 2 ( e) � 1, (31) and (32) as well as the controller yields

which is negative definite and guarantees that the trajectories of the state vector ( e, Z) T converge asymptotically to the origin.

B. Kinematic controller

We proceed to design the controller to head the quad rotor velocity Vr E D in such a way that the angular error is rendered to zero, which means that the desired path has been achieved. Unlike classical kinematic models where the heading rate w p is interpreted as the actual control input, in the present case, however, the control input is the torque T 1j;.

Hence, the kinematic subsystem is rewritten as The subsystem (38-39) is stabilized through a hierarchical control scheme (two-time scale), considering the path devi ation kinematics as the slow-time scale (outer-loop) and the heading dynamics as the fast time-scale (inner-loop). Since the quad-rotor kinematics is driven via the angle 'ijJ, which is not a control but a state variable of the heading dynamics.

Then, the heading is controlled by T1j; to track the stabilizing


Marwa Ben Jabra

Ramzi Guetari

Aladine Chetouani

Hedi Tabia

Nawres Khlifa

Facial expression recognition using the bilinear pooling

Keywords: Facial Expression Recognition, Image Classification, Deep Learning, Bilinear Pooling, Bilinear-CNN

INTRODUCTION

The automatic recognition of humans' facial expressions is an emerging field in scientific research. The latter allows new applications in other areas such as human-computer interaction, security, medicine, etc. The goal of a facial expressions recognition system is to identify simple or combined universal emotions by analyzing the characteristics of facial deformities of the mouth, eyes and eyebrows. Reliable recognition of facial expressions under natural conditions is therefore necessary, but remains an area where many problems are still unresolved (variations in lighting, position of the head, occlusions, etc).Artificial Intelligence has made a quantum leap in recent years, with the advent of neural networks and machine learning. These techniques, combined with computers, sensors and cameras, have produced increasingly intelligent and autonomous systems. They have also made it possible to significantly change the artificial vision that is at the base of the recognition of individuals and their facial expressions. In the literature, several approaches have been proposed. Recently a special attention has been given to the methods based on the deep learning and particularly CNNs. The latter showed a real potential to overcome certain difficulties in the process of recognition and have thus significantly increased the classification rate.

Several works are devoted to this field but differ considerably by the CNN architectures adopted, the learning and test protocols applied, the parameters used (filter size, number of filters, etc.), etc. Indeed, the results that we studied in the literature and experimented with different architectures differ considerably from one model to another, depending on the choice of parameters. In addition, the same model can produce different results on the same set of data with two different initializations.

The fact that two CNNs give different results or that the same CNN with different initialization can give different results has spurred the idea of using the bilinear pooling (Liu et,2014) that consists in using two CNNs in parallel as described in § C. We believe that bilinear CNNs drastically improve the recognition performance comparing to previous works. This is what we are trying to demonstrate in this paper.

The contributions we bring through this study can be summarized as follows: first of all, we propose the use of a pre-trained model (such as VGG16, Inception or ResNet50) with a transfer learning to classify the basic emotions. We then propose our own shallow CNN architecture built from scratch for the same purpose. This architecture has fewer layers than predefined architectures and a smaller number of parameters. Finally, we propose the use of the bilinear pooling to improve the results.

The rest of this paper is organized into four sections. In Section 2, we present the related work on the recognition of facial emotions with the CNNs and the bilinear CNN models, as well as the relevant literature on the combination of features. Then, In Section 3, we describe our contribution and our proposed deep neural network architectures for emotion recognition. In Section 4, we present our experimental results. The last section is devoted to the conclusion.

RELATED WORK

The recognition of facial emotions is based on facial landmarks. The most significant facial references that have a real impact on the facial analysis and recognition are the corners of the eyes, the eyebrows, the lobes of the ears, the tip of the nose, nostrils, chin, mouth, etc. While humans are naturally able to understand these emotions, facial landmarks identification is a challenging task in automatic emotions recognition. In fact, the performance of an automated approach depends on the selection of a discriminate feature set and the use of an efficient learning technique. Recent researches have shown the potential and effectiveness of the use of CNNs. Unlike other ANNs that require a pre-processing phase to set the learning features, CNNs are deep artificial neural networks in which, the convolution blocks automatically select the most relevant filters to extract relevant features.

Behzad Hasani et al. [START_REF] Behzad | Facial Expression Recognition Using Enhanced Deep 3D Convolutional Neural Networks[END_REF] proposed a 3D CNN architecture named 3D Inecption-ResNet (3DIR) for the extraction of spatial relationships in facial images and the temporal relationships between different video frames. The approach extracts 66 landmarks from the face using a regression method of local binary characteristics [START_REF] Dhananjay | Facial Landmark Localization -A Literature Survey[END_REF]. These landmarks are used as inputs to the 3DIR model during the training phase. Toevaluate the proposed method they used four facial expression databases, which are MMI [START_REF] Pantic | Webbased database for facial expression analysis[END_REF], CK+ [START_REF] Lucey | The extended cohn-kanade dataset (ck+): A complete dataset for action unit and emotionspecified expression[END_REF], GEMEP-FERA [START_REF] Banziger | Introducing the geneva multimodal emotion portrayal (gemep) corpus. Blueprint for affective computing: A sourcebook[END_REF], and DISFA [START_REF] Mavadati | Disfa: A spontaneous facial action intensity database[END_REF]. They obtained a recognition rate of 67.52% on CK+ dataset.

Zhan Wu et al [START_REF] Zhan | NIRExpNet: Three-Stream 3D Convolutional Neural Network for Near Infrared Facial Expression Recognition[END_REF] proposed a threestream 3D CNN that automatically extracts both spatial and temporal characteristics. The purpose is to fuse local and global facial expression features by using different methods that could produce different recognition rates. The concatenation method has given the best recognition rate of 78.42%.

Octavio Arriaga et al. (Arriaga 2017) proposed two real-time CNN models to simultaneously perform the face detection, the gender recognition and the emotion classification tasks. In the first model, they removed the fully connected layers and used a Global Average Pooling, which reduces each feature map to a scalar value. The second model is inspired by the Xception [START_REF] Chollet | Xception: Deep Learning with Depthwise Separable Convolutions[END_REF] architecture in which fully connected layers are eliminated and residual modules [START_REF] Kaiming | Deep residual learning for image recognition[END_REF]) have been included. They also proposed to use depth-wise separable convolutions that separate the processes of features extraction to depth-wise convolution layers and point-wise convolution layers to separate the spatial cross-correlations from the channel crosscorrelations. The accuracy rate for emotion recognition and classification reached 66% using the FER-2013 dataset.

Liu et al. (Liu 2014) proposed a CNN model named 3DCNN-DAP for the emotion recognition, based on Deformable Action Parts (DAP). This work focuses on the deformable facial action part, which is the most important facial information when analyzing emotions and dynamically changes according to the facial expressions. They proposed a technique for learning deformable parts of action and for detecting specific parts of the action of the face under structured spatial constraints. To evaluate the proposed method, they used CK+ and MMI posed expression datasets and, the FERA spontaneous dataset. The best accuracy rate they have obtained is 87.5% on CK+ dataset.

In their paper, Khorrami et al. [START_REF] Khorrami | Do Deep Neural Networks Learn Facial Action Units When Doing Expression Recognition[END_REF]) proposed a CNN architecture for the selection of the most relevant parts of the human face allowing the recognition of facial expressions. They built a zerobias CNN network, which consists on ignoring the biases of the convolutional layers in order to reduce the number of parameters. To evaluate the performance of this method, they used the CK + and Toronto Face Dataset (TFD) datasets. They obtained an accuracy rate of 81.8% with the CK + dataset and an accuracy rate of 79.4% with the TFD dataset.

Li et al. [START_REF] Li | Patch-Gated CNN for Occlusion-aware Facial Expression Recognition[END_REF] proposed a Patch-Gated CNN (PG-CNN) model that breaks down the face into facial regions. They extracted an intermediate features map of the facial image using VGGnet model. After, this intermediate features map is decomposed by the PG-CNN into 24 sub-featuremaps corresponding to 24 local patches. They then calculated the weight of each patch to perceive facial occlusions. Finally, the representation of the occluded face is obtained by the concatenation of the weighted local features. They achieved an accuracy of 80.28% using the CK + database.

In the paper of Xie and Hu [START_REF] Xie | Facial expression recognition using hierarchical features with deep comprehensive multipatches aggregation convolutional neural networks[END_REF], the Deep Comprehensive Multi-patches Aggregation CNN (DCMA-CNN) for the recognition of facial expression has been introduced. This architecture is composed of two branches of a CNN. One branch extracts local features depicting the expression details of image patches, while the other extracts global features characterizing the high-level semantic information of an expression derived from the entire expression image. The local and global features are fused for the classification task. To evaluate the model they used CK+ and JAFFE facial expression datasets, and they obtained a rate accuracy of 88.67%.

Zhou et al. [START_REF] Zhou | Fine-Grained Facial Expression Analysis Using Dimensional Emotion Model[END_REF], proposed a method that uses bilinear-CNN models. They transformed facial expression analysis from a classification problem into a regression one to predict the facial expression updates based on local appearance features. They explored two different deep CNN architectures, VGG16 and ResNet50. To improve these methods by a bilinear pooling, they modified them accordingly to serve the regression task. Specifically, they removed all the layers after the global pooling layer that aggregated all features spatially; these layers include the last convolution layer and the loss layer. Over the global pooling layer, they stacked a new convolution layer with randomly initialized weights to output final embedding vectors for input images. They evaluated their models on the FER-2013 dataset and they obtained a rate accuracy of 81,79% using global pooling VGG16 and 83,78% using bilinear pooling VGG16. The proposed method recognize the basic emotion and a full range of other emotions.

Zhang et al. [START_REF] Zhang | Deep Fusion: An Attention Guided Factorized Bilinear Pooling for Audio-video Emotion Recognition[END_REF] proposed a Bilinear CNN model based on factorized bilinear pooling (FBP) that aims to use the audio and video features for the emotion recognition. In this proposed method, they extracted the emotion features from video and audio using two different parallel streams of CNN. Then, they fused the outputs feature maps using factorized bilinear pooling. Finally, they classify the emotion using a fully connected layer followed by a Softmax function. To validate their method, they used the audio-video AFEW dataset and they obtained a rate accuracy of 62,48%.

PROPOSED METHOD

The main contribution of this work is the use of bilinear pooling by combining a standard CNN (with generally performs well, but suffers from a huge number of parameters and are time consuming) and a shallow CNN built from scratch, which has a reduced number of parameters, and an acceptable recognition rate. Bilinear pooling is considered as a Second-order aggregation of CNN activations, which can provide improvements over classical representations using first-order aggregation (e.g., sum or max). Actually, the Bilinear CNN (B-CNN) is considered as a stateof-the-art network architecture for texture and finegrained recognitions. A bilinear CNN model consists of two parallel CNNs, each of them extracts the features map from an input image. Features maps produced by each of the two CNNs are multiplied using the outer product generating second-order characteristics. The latter are then pooled to form high dimensional bilinear features to obtain an image descriptor. Bilinear features maps are normalized using sqrt (square root) and L2 normalization. Finally, fully connected layers and a Softmax function are used for the classification task. Such bilinear architectures have been proposed for the recognition and classification of facial emotions. Bilinear CNN models have proven effective for fine recognition, scene categorization, texture recognition, and visual question-and-answer tasks, among others. They are able to distinguish the subtle differences between cars, birds and planes [START_REF] Lin | Bilinear Convolutional Neural Networks for Fine-grained Visual Recognition[END_REF] (Lin et al, 2015). In the following, we present the transfer learning process, the architecture of the shallow CNN as well as the bilinear architecture.

Transfer Learning for Emotion Recognition

First, we proposed to explore standard CNN architectures such as (Densenet, Mobilenet V1 & V2, InceptionResnetV2, NasnetMobile, Xception, ResNet50, VGG16 and VGG19), which are widely studied and exploited in practical applications. To train these architectures to the classification of the basic facial expression, we combined the transfer learning and the fine tuning. The transfer learning is a machine learning technique that consists of using a model trained in a particular task to perform another one. The fine tuning consists in using an already trained model and training it on a specific dataset. This technique requires that the last fully connected layer be replaced by another one adapted to the new task to be performed (classifying the seven basic facial expressions in our case). These deep models have performed well, but they use a huge number of parameters for the training.

Shallow CNN Architecture

To obtain a good performance with small number of parameters, we proposed to introduce a new CNN architecture with a small number of layers. After several experimentations, we built a shallow CNN with 7 convolution layers, 3 max-pooling layers and a single fully connected layer. The table 1 summarizes the architecture of the proposed shallow CNN. To train this architecture to the classification of the basic facial expression, we used transfer learning and the fine tuning. Our architecture has been trained on facial expressions datasets. The classification of the facial expressions with the shallow architecture has given the accuracy rate of 67.78% (on the CK+ dataset), which is quite good considering the depth of the CNN and its number of parameters.

Bilinear Pooling for Emotion Recognition

A bilinear CNN model consists of two parallel CNNs, each of them extracting the features map from an input image. Features maps produced by each of the two CNNs are multiplied using the outer product producing second-order characteristics. The latter are then pooled to form high-dimensional bilinear features to obtain an image descriptor. Bilinear features maps are normalized using sqrt (square root) and L2 normalization. Finally, fully connected layers and a Softmax function are used for the classification task. Such bilinear architectures have been proposed for the recognition and classification of facial emotions. Bilinear CNN models have proven effective for fine recognition, scene categorization, texture recognition, and visual question-and-answer tasks, among others. They are able to distinguish the subtle differences between cars, birds and planes (Lin 2017) (Lin 2015).

The bilinear architecture nevertheless poses some performance problems in terms of temporal complexity. The learning phase of both networks is time consuming. Even the recognition phase can often be time consuming. This is due to the depth of the architectures used and the number of parameters required for each branch of the bilinear architecture. Our goal is to reduce the processing time of a bilinear architecture while keeping a very good accuracy rate. The combination of a standard architecture with our shallow architecture seemed to us a good compromise that needed to be tested.

We tested a number of bilinear configurations. We combined two instances of standard architectures as well as combined with our shallow CNN. All the CNNs have been modified to become a bilinear architecture. All layers after the last convolution layer have been removed, namely the last pooling layer, the fully connected layer, and the loss layer. The weights of each CNN have also been fixed. The two CNNs were used in parallel with a different weight initialization for each of them in order to extract the feature maps from the same input images. These feature maps were merged using the Khroneker product [START_REF] Paumard | Image Reassembly Combining Deep Learning and Shortest Path Problem[END_REF]. Two normalization layers follow the features extraction (sqrt and L2 have experimented with more or less the same performance). Finally, a fully connected layer with 7 outputs followed by a Softmax function was used for the final classification. The bilinear architecture that combines VGG16 and our shallow CNN is presented in Fig. 1.

RESULTS AND DISCUSSION

In this section, we present the facial expression databases used in our experiments and we detail the experimental results.

Results

We used the same CNN standard architecture with different initializations and we combined them to form the bilinear CNN model. The accuracy rate using the bilinear architecture has been improved compared the initial results. The bilinear VGG16 model gives the best accuracy rate but it still uses a huge number of parameters. A bilinear architecture using two instances of the shallow CNN has also been tested. It resulted in the third best accuracy while using a reduced number of parameters. Finally, we implemented a bilinear architecture by combining VGG16 with the shallow CNN to see the impact on the accuracy rate while reducing the number of parameters and then reduce the processing time.

Table 2 and table 3 summarize the results obtained using standard configurations and bilinear ones on respectively CK+ and FEI dataset. For the CK+ dataset, the fine-tuned VGG16 achieved the best accuracy rate, while the shallow CNN has achieved the third performance with an accuracy rate of 67,78%. In a bilinear configuration, VGG16 has also achieved the best performance with an accuracy rate of 84.30 %, and the shallow CNN has realized an acceptable performance with an accuracy rate 78,65%.

The combination of VGG16 and the shallow CNN achieved the best performance of all standard and bilinear configurations tested. It achieved the accuracy rate of 86,98%.For the, FEI Dataset, the VGG16 in a standard configuration reached the accuracy rate of 80.89%, which is the best performance of all the tested architectures. Our shallow CNN has given the third performance with an accuracy rate of 67,20%. The VGG16 in a bilinear configuration has given the best performance of all the standard architectures used in the same configuration. It, indeed reached the accuracy rate of 83.52%. Our shallow model used in a bilinear configuration has given the third accuracy rate that is 77,15%. Finally, the bilinear configuration in which we combined the VGG16 model with our shallow model gives the best performance of all the bilinear configurations as well as he standard configuration with an accuracy rate of 85,35%.

Comparison with the State of the Art

In In Table 5, we compare our proposed CNN models with (PCA (Wang 2013), LDA methods (Thomaz 2010) using FEI dataset. We notice that the handcrafted methods implemented FEI dataset give an accuracy rate better than the proposed model.

For the fine-grained facial expression analysis, we experimented various CNN architectures such as VGG16, ResNet50, Mobilnet, as well as our own shallow model. We evaluated these architectures on a standard CK + and FEI databases. The results have shown that deep CNN models achieve very promising results. Moreover, the bilinear CNN models have significantly outperformed their respective counterparts.

CONCLUSIONS

Advances in modern techniques in artificial intelligence have led to considerable advances in the field of facial recognition. Nevertheless, there are still a number of problems to solve in order to achieve performances close to those of humans. Our contribution is, in our opinion, a contribution to advancing research in this area in the right direction. Normally, the use of CNNs imposes the introduction of a large number of parameters with more or less efficient classification rates. The most powerful architectures often use the largest number of parameters. We tested these architectures individually and combined using the bilinear pooling. We found that a bilinear architecture improves performance, but the number of parameters remains too large.

To measure the impact of the number of parameters in the classification performance, we developed our own shallow architecture which, despite the small number of parameters, still gave good results when it's used individually.

We then used our shallow model with the standard VGG16 architecture that gave us the best performance. We evaluated our proposed bilinear model on the CK + and FEI datasets. We found that this model allowed us to obtain a good precision rate, namely 86.98% using the CK+ dataset, and 85.35% on the FEI dataset while using a reduced number of parameters. We can therefore conclude that the bilinear model CNN has clearly exceeded the performance of simple deep CNN models.

We compared the results with other deep models and handcrafted methods. We found that we performed better than some methods. We also found that methods that use local and global characteristics exceeded our method in terms of performance.

In perspective, we will propose a CNN model that extracts local and global characteristics and merges them with bilinear pooling to improve the performance of our method.

Figure 1 :

 1 Figure 1: The proposed Bilinear CNN architecture.

Table 1 :

 1 Architecture of the shallow CNN.

	Layer	Output shape	Params. count
	Input Layer	224, 224, 3	0
	Convolution	224, 224, 64	1792
	Convolution	224, 224, 64	36928
	Max-pooling	112, 112, 64	0
	Convolution	112, 112, 128	73856
	Convolution	112, 112, 128	147584
	Max-pooling	56, 56, 128	0
	Convolution	56, 56, 256	295168
	Convolution	56, 56, 256	590080
	Convolution	56, 56, 256	590080

Table 2 :

 2 Results obtained using standard configurations and bilinear ones on CK+ dataset.

			Standard CNNs			Bilinear pooling
	ARCHITECTURES	Training Accuracy	Test accuracy	Number of parameters	Training Accuracy	Test accuracy	Number of parameters
	DenseNet (1) DenseNet (2)	53,10% 52,68 %	35,96% 33,40%	8,062,504	91,08%	34,09%	7,040,538
	MobileNet-V1 (1) MobileNet-V1 (2)	53,75% 56,06 %	31,09% 28,46%	4,253,864	96,53%	37,38%	14,568,903
	MobileNet V2 (1) MobileNet V2 (2)	60,51 % 61,32 %	29,59% 31,87%	3,538,984	100,00%	42,82%	13,721,671
	NASNetMobile (1) NASNetMobile (2)	62.10% 63,30 %	30,90% 34,09%	5,326,716	75,43%	41,18%	8,526,964
	InceptionResNetV2 (1) InceptionResNetV2 (2)	61,5 7% 59,77 %	36,13% 33,41%	55,873,736	82,62%	39,20%	49,758,951
	ResNet50 (1) ResNet50 (2)	58,38% 59,03 %	40,07% 42,98%	25,636,712	86,90%	41,50%	44,363,911
	Xception (1) Xception (2)	51,68% 53,99 %	42,31% 44,01%	22,910,480	77,48%	36,77%	50,221,615
	VGG16 (1) VGG16 (2)	62,21% 65,89 %	69,16% 70,60 %	138,357,544	100,00%	84,30%	16,549,703
	VGG19 (1) VGG19 (2)	63,91% 64,88 %	68,91% 69,33%	143,667,240	100,00%	83,12%	21,859,399
	Shallow (1) Shallow (2)	65,21% 62,97 %	67,78 % 65,77%	3,140,423	100,00%	78,65%	3,929,735
	VGG16 (2) Shallow (1)	65,89 % 65,21%	70,60% 67,78%	143,667,240 3,140,423	100,00%	86,98%	9,470,279

Table 3 :

 3 Results obtained using standard configurations and bilinear ones on FEI dataset.

	ARCHITECTURES	Training Accuracy	Validation accuracy	Number of parameters	Training Accuracy	Validation accuracy	Number of parameters
	DenseNet (1) DenseNet (2)	60.34% 59,63 %	25,32 % 23,99 %	8,062,504	80,52%	31,09%	7,040,538
	MobileNet-V1 (1) MobileNet-V1 (2)	60.62% 61,08 %	29,63 % 30,10 %	4,253,864	68,55%	31,62%	14,568,903
	MobileNet V2 (1) MobileNet V2 (2)	64,75 % 65,03 %	32,23 % 34,73 %	3,538,984	75,65%	42,33%	13,721,671
	NASNetMobile (1) NASNetMobile (2)	54.03% 57,88 %	22,91 % 25,83 %	5,326,716	65,43%	32,53%	8,526,964
	InceptionResNetV2 (1) InceptionResNetV2 (2)	61,5 7% 60,55 %	36,13% 34,11 %	55,873,736	81,43%	45,20%	49,758,951
	ResNet50 (1) ResNet50 (2)	66,85% 63,36 %	50,49% 48,77 %	25,636,712	86,90%	51,31%	44,363,911
	Xception (1) Xception (2)	51,68% 49,65 %	32,31% 36,98 %	22,910,480	57,58%	45,46%	50,221,615
	VGG16 (1) VGG16 (2)	70,09% 69,83 %	80,89 % 79,89 %	138,357,544	100,00%	83,52%	16,549,703
	VGG19 (1) VGG19 (2)	68,39% 66,97 %	78,21% 76,88%	143,667,240	100,00%	82,77%	21,859,399
	Shallow (1) Shallow (2)	68,20% 65,37 %	67,20 % 64,20 %	3,140,423	97,67%	77,15%	3,929,735
	VGG16 (2) Shallow (1)	70,09 % 67,21%	80,89 % 57,33 %	143,667,240 3,140,423	100,00%	85,35%	9,470,279
	4.1 Databases			it is composed of 327 video sequences. We extracted
				the last three images from each sequence of the CK+
	In our experiments, we used CK+ and FEI facial	video dataset to build a new image dataset with 981
	expression datasets, which consist of seven basic	mainly grey images. The FEI facial expression
	expressions (anger, disgust, fear, happiness, sadness	dataset has been created in the Artificial Intelligence
	and surprise). The CK+ dataset (Cohn-Kanade	Laboratory of FEI Brazil and it is composed of 252
	extended dataset) is the extension of the CK one and	images.

Table 4

 4

	, we compare the performance of the
	proposed Bilinear model on CK+ dataset, with
	different other methods, including Hand-crafted-
	based methods (TMS (Lin 2017), CDA (Jain 2011),

Table 4 :

 4 Results of hand-crafted and CNN models on (CK+) Dataset.

		Method	Validation accuracy
		TMS (Lin 2017)	91.80%
	Hand-crafted	CDA (Jain 2011)	85.00%
	features	MSR (Jain 2011)	91.40%
		ITBN (Rifai 2012)	86.30%
		3DCNN-DAP (Kaiming 2016)	87.90%
	CNN	Zero-bias CNN (Liu 2014b)	81.80%
	Architecture	PG-CNN (Khorrami 2015)	80.28%
		LFCNN (Zhao 2011)	88.67%
		VGG16	70.60%
	Our CNN models	Shallow CNN Bilinear VGG16 Bilinear Shallow CNN	67.78% 84.30% 78.65%
		Proposed Bilinear model	86.98%

Table 5 :

 5 Results of hand-crafted and CNN models on FEI dataset.

		Architectures	CNN	Proposed Bilinear CNN
	CNN Architecture	VGG16 Shallow CNN	70,0% 67,2%	85,35%
	Handcrafted	PCA (Wang 2013)		94,00%
	features	LDA (Thomaz 2010)		92,20%