Vincent Drouard

Radu Horaud

Antoine Deleforge

Silèye Ba

Georgios Evangelidis

Robust Head-Pose Estimation Based on Partially-Latent Mixture of Linear Regressions

Keywords: Head pose, face detection, mixture of linear regressions, manifold learning, expectation-maximization

Head-pose estimation has many applications, such as social event analysis, human-robot and human-computer interaction, driving assistance, and so forth. Head-pose estimation is challenging because it must cope with changing illumination conditions, variabilities in face orientation and in appearance, partial occlusions of facial landmarks, as well as bounding-boxto-face alignment errors. We propose tu use a mixture of linear regressions with partially-latent output. This regression method learns to map high-dimensional feature vectors (extracted from bounding boxes of faces) onto the joint space of head-pose angles and bounding-box shifts, such that they are robustly predicted in the presence of unobservable phenomena. We describe in detail the mapping method that combines the merits of unsupervised manifold learning techniques and of mixtures of regressions. We validate our method with three publicly available datasets and we thoroughly benchmark four variants of the proposed algorithm with several state-of-the-art head-pose estimation methods.

I. INTRODUCTION

Head pose is an important visual cue in many scenarios such as social-event analysis [START_REF] Sabanovic | Robots in the wild: observing human-robot social interaction outside the lab[END_REF], human-robot interaction (HRI) [START_REF] Goodrich | Human-robot interaction: A survey[END_REF] or driver-assistance systems [START_REF] Murphy-Chutorian | Head pose estimation for driver assistance systems: A robust algorithm and experimental evaluation[END_REF] to name a few. For example, in social-event analysis, 3D head pose information drastically helps to determine the interaction between people and to extract the visual focus of attention [START_REF] Murphy-Chutorian | Head pose estimation in computer vision: A survey[END_REF]. The pose is typically expressed by three angles (pitch, yaw, roll) that describe the orientation with respect to a head-centered frame. The estimation of the pose parameters is challenging for many reasons. Algorithms for head-pose estimation must be invariant to changing illumination conditions, to the background scene, to partial occlusions, and to inter-person and intra-person variabilities. In most application scenarios, faces have small support area, i.e. bounding boxes, typically of the order of 100×100 pixels. Even if the face bounding box is properly detected, one has to extract the pose angles from low-resolution data.

Recent advances in computer vision have shown the relevance of representing an image patch with a feature vector, e.g. SIFT [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], HOG [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF], SURF [START_REF] Bay | Speeded-up robust features (surf)[END_REF], or one of their variants. The rationale of representing faces in a high-dimensional feature This work has been funded by the European Research Council through the ERC Advanced Grant VHIA #340113. R. Horaud acknowledges support from a XEROX University Affairs Committee (UAC) grant (2015-2017).

space is that the latter supposedly embeds a low-dimensional manifold parameterized by the pose parameters, or the headpose manifold, e.g. [START_REF] Osadchy | Synergistic face detection and pose estimation with energy-based models[END_REF]- [START_REF] Deleforge | High-dimensional regression with gaussian mixtures and partially-latent response variables[END_REF]. Hence, several attempts were carried out in order to cast the problem at hand into various frameworks, such as manifold learning (unsupervised) [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF], regression [START_REF] Deleforge | High-dimensional regression with gaussian mixtures and partially-latent response variables[END_REF], [START_REF] Drouard | Head pose estimation via probabilistic high-dimensional regression[END_REF], convolutional neural networks [START_REF] Osadchy | Synergistic face detection and pose estimation with energy-based models[END_REF] (supervised), or dimensionality reduction followed by regression [START_REF] Foytik | A two-layer framework for piecewise linear manifold-based head pose estimation[END_REF], to cite just a few.

While the papers just cited yield interesting and promising results, there are several major issues associated with representing faces with high-dimensional feature vectors, issues that have not been properly addressed. Indeed and as already mentioned, these vectors contain many underlying phenomena other than pose, e.g. illumination, appearance, shape, background, clutter, etc. Hence, one major challenge is to be able to remove these other pieces of information and to retain only head-pose information. Another drawback is that head pose relies on face detection, a process that amounts to finding a bounding box that contains the face and which is invariant to face orientation.

Take for example the case of finding pose parameters using linear or non-linear manifold learning followed by regression. This is usually justified by the fact that high-dimensionalto-low-dimensional (high-to-low) regression has to estimate a very large number of parameters, typically of the order of D 2 where D is the dimension of the feature space. This in turn requires a huge training dataset. Moreover, this sequential way of doing presents the risk to map the input onto an intermediate low-dimensional space that does not necessarily contain the information needed for the finally desired output -head pose. Finally, the estimation of the pose angles in two distinct steps cannot be conveniently expressed in a single optimization formulation. Supervised pose estimation also requires an annotated dataset of faces with their bounding-box locations and the corresponding pose parameters. The large majority of existing methods relies on manually annotated faces for training. In particular this ensures good bounding-box-to-face alignment, i.e. the face center coincides with the bounding-box center and the bounding box does not contain too many background pixels. This alignment requirement has, however, an important drawback, namely that bounding boxes found by face detection techniques do not correspond to the annotated ones. This means that feature vectors obtained with a face detector do not Fig. 1. Pipeline of the proposed supervised head-pose estimation method. Top: the parameters of a mixture of linear regressions are learnt from faces annotated with their poses (left). The result of this learning is a simultaneous partitioning of both the high-dimensional input (high-dimensional feature vectors shown in the middle) and low-dimensional output (two-dimensional parameter space shown on the right), such that each region in this partition corresponds to an affine mapping between the input and the output. Moreover, the output is modeled by a Gaussian mixture and each region corresponds to a mixture component. This yields a predictive distribution that can then be used to predict an output from a test input. Bottom: A face detector is used to localize a bounding box (left, shown in red) from which a HOG descriptor, namely a high-dimensional feature vector, is extracted. Using the predictive distribution just mentioned, it is then possible to estimate the head-pose parameters (yaw and pitch in this example). Additionally, it is also possible to refine the bounding-box location such that the latter is optimally aligned with the face (right, shown in green).

necessarily belong to the space of feature vectors obtained by manual annotation and used for training. Hence, there is a risk that the head pose, predicted from feature vectors associated with automatic face detection, is incorrectly estimated.

In this paper we propose to learn with both head-pose parameters and bounding-box-to-face alignments, such that, at runtime both the head-pose angles and bounding-box shifts are predicted. This ensures that the bounding-box-to-face alignments predicted with our method are similar with those used for training. Consequently, background variations have minimal influence on the observed feature vector from which the pose is being predicted. This prevents from pose-estimation errors due to discrepancies that occur between manually annotated faces (used for training) and automatically detected faces (used for prediction). We adopt a high-dimensional description of faces based on HOG features, hence we need to solve a high-to-low regression problem. Rather then performing dimensionality reduction and regression in two distinct steps, we propose to use a generative model that unifies these two steps in a single one.

More precisely, we propose to adopt the partially-latent mixture of linear regression model of [START_REF] Deleforge | High-dimensional regression with gaussian mixtures and partially-latent response variables[END_REF]. When applied to head pose estimation from feature vectors, this method has two advantages: (i) it solves high-dimensional to low-dimensional regression in a single process, without the need of a dimensionality reduction pre-processing step, and (ii) the method incorporates latent variable augmentation. More precisely, the output variables are only partially observed such that the regression function can be trained with partially observed outputs. The method adopts an inverse regression strategy [START_REF] Li | Sliced inverse regression for dimension reduction[END_REF]: the roles of the input and output variables are interchanged, such that high-to-low regression is replaced with low-to-high regression, followed by high-to-low prediction. As illustrated by our experiments, e.g. Fig. 1, we found that the inclusion of a partially-latent output is particularly relevant whenever the high-dimensional input is corrupted by irrelevant information that cannot be easily annotated, e.g. facial expressions, changes in illumination conditions, changes in background, etc.

Another contribution is a thorough experimental evaluation and validation of the merits of [START_REF] Deleforge | High-dimensional regression with gaussian mixtures and partially-latent response variables[END_REF] when the task is to recover low-dimensional parameterizations of manifolds embedded in high-dimensional feature vector representations of image patches. We note that this task is ubiquitous in computer vision and image processing, hence a method that unifies dimensionality reduction and regression in a principled manner is of primary importance and of great utility for practitioners. Moreover, with respect to [START_REF] Deleforge | High-dimensional regression with gaussian mixtures and partially-latent response variables[END_REF] we provide more detailed mathematical derivations as well as a unified view of the algorithms used for training, such that the method can be readily reproduced by others. 1The remainder of the paper is organized as follows. Section II discusses related work. The regression method is described in Section III and the algorithmic details can be found in Appendix A. Section IV describes implementation details. Experimental results are presented and discussed in Section V. Finally, Section VI draws some conclusions for future work.

II. RELATED WORK

Head-pose estimation has been very well investigated for the past decades and several powerful methods have been developed. The head-pose estimation literature was surveyed a few years ago [START_REF] Murphy-Chutorian | Head pose estimation in computer vision: A survey[END_REF], however this survey does not include methods based on depth data as these papers were published after 2009. For convenience we grouped the head-pose methods into four categories: (i) methods that are based on depth images, (ii) methods based on manifold learning, (iii) methods based on regression, and (iv) methods that combine pose with face detection and localization.

The recent advent of depth cameras enabled fast development of depth-based head-pose methods. Depth data allow to overcome some of the drawbacks of RGB data, such as illumination problems and facial landmarks detection, which is more reliable. One of the first methods that uses depth data is [START_REF] Seemann | Head pose estimation using stereo vision for human-robot interaction[END_REF], were a depth map of the head region is combined with a color histogram and used to train a neural network. Random forest regression is proposed in [START_REF] Fanelli | Random forests for real time 3d face analysis[END_REF] to estimate both head pose and facial landmarks. Random forest regression is also used in [START_REF] Peng | From circle to 3-sphere: Head pose estimation by instance parameterization[END_REF] where RGB SIFT descriptors are combined with 3D HOG descriptors. More recently [START_REF] Ghiass | Highly accurate and fully automatic head pose estimation from a low quality consumer-level rgb-d sensor[END_REF] released a new method combining RGB and depth images to infer head pose. The RGB image is used to find facial landmarks, the 3D positions of the landmarks are used to fit a plane on the face that will be used to extract the 3D points that belong to the face. Using the face 3D point clouds, a morphable model of a face is map on it using optimization methods to estimate head orientation parameters. A general observation about these methods is that depth information is merely used to disambiguate photometric data and that depth data cannot be used alone for head pose.

Several authors proposed to use manifold learning, namely finding a low-dimensional output space of head poses from a high-dimensional input space of feature vectors. Nevertheless, the output variables do not necessarily correspond to the pose angles, so one has to farther learn in a supervised way the mapping between the manifold-learning output and the desired space spanned by the pose parameters. This has been achieved in various ways, e.g. [START_REF] Foytik | A two-layer framework for piecewise linear manifold-based head pose estimation[END_REF], [START_REF] Srinivasan | Head pose estimation using view based eigenspaces[END_REF]- [START_REF] Sundararajan | Head pose estimation in the wild using approximate view manifolds[END_REF]. As already mentioned, these two step methods suffer from the fact that unsupervised manifold-learning techniques do not guarantee that the predicted output space contains the information needed for head pose.

Among the regression methods used for head pose are Gaussian process regression (GPR) [START_REF] Marin-Jimenez | Detecting people looking at each other in videos[END_REF], support vector regression (SVR) [START_REF] Murphy-Chutorian | Head pose estimation for driver assistance systems: A robust algorithm and experimental evaluation[END_REF], partial least squares (PLS) [START_REF] Sharma | Bypassing synthesis: Pls for face recognition with pose, low-resolution and sketch[END_REF] and and kernel PLS [START_REF] Haj | On partial least squares in head pose estimation: How to simultaneously deal with misalignment[END_REF]. Both [START_REF] Marin-Jimenez | Detecting people looking at each other in videos[END_REF] and [START_REF] Murphy-Chutorian | Head pose estimation for driver assistance systems: A robust algorithm and experimental evaluation[END_REF] estimate the pose angles independently, so several regression functions must be learned, one for each angle, hence correlations between these parameters cannot be taken into account. Another drawback of all kernel methods is that they require the design of a kernel function with its hyper-parameters, which must be either manually selected or properly estimated using non-convex optimization techniques.

PLS and kernel PLS proceed in two steps. First, both the input and the output are projected onto low-dimensional latent subspaces by maximizing the covariance between the projected input and the projected output. Second, a linear regression between these two latent subspaces is estimated. The performance of PLS is subject to the relationship between the covariance matrices of input and output variables and to the eigen structure of the covariance of the input variable [START_REF] Naik | Partial least squares estimator for single-index models[END_REF]. The advantage of the proposed method is that it estimates a mixture of linear regressions directly from the input and output variables.

The methods described so far use manually annotated images for training and face detectors for testing. As already discussed, this could lead to pose estimation errors because a test feature vector may not lie in the subspace spanned by the feature vectors used for training. One way to deal with this problem is to combine face detection and pose estimation in a single task. For example, [START_REF] Zhu | Face detection, pose estimation, and landmark localization in the wild[END_REF] considers face templates for all possible poses which are then fed into a cascaded classifier.

Convolutional neural network (CNN) architectures were also proposed in the recent past [START_REF] Osadchy | Synergistic face detection and pose estimation with energy-based models[END_REF], [START_REF] Ahn | Real-time head orientation from a monocular camera using deep neural network[END_REF]. [START_REF] Osadchy | Synergistic face detection and pose estimation with energy-based models[END_REF] considers a fixed image sub-window at all locations and scales. The network consists of 64, 000 weights and kernel coefficients that need to be estimated, and both face and non-face samples must be considered. Altogether, training the network with 52, 000 positives and 52, 000 negatives involves non-linear optimization and takes 26 hours on a 2GHz Pentium 4. [START_REF] Ahn | Real-time head orientation from a monocular camera using deep neural network[END_REF] proposed a CNN architecture composed of four convolutional layers with max-pooling on the first two layers; the activation function is the hyperbolic tangent which yields good convergence during the training phase. Small input RGB images (32 × 32 pixels) and small filters (5 × 5 pixels) are used in order to overcome the limitation of the training dataset. The network is trained using 13, 500 face patches extracted from the dataset. More recently, [START_REF] Liu | 3D head pose estimation with convolutional neural network trained on synthetic images[END_REF] proposed to simulate a dataset of head poses in order to train a CNN. Then they use the trained network to estimate head pose from real color images using the BIWI dataset [START_REF] Fanelli | Random forests for real time 3d face analysis[END_REF]. They show that when trained using the BIWI dataset the CNN approach yields results similar to [START_REF] Drouard | Head pose estimation via probabilistic high-dimensional regression[END_REF] and that the accuracy is improved, by a factor of 2, when a large set of simulated images are used for training.

The problem of bounding-box-to-face miss-alignment was discussed and addressed in [START_REF] Haj | On partial least squares in head pose estimation: How to simultaneously deal with misalignment[END_REF]. First, kernel PLS is trained on manually extracted bounding boxes and associated pose parameters, and second, head-pose prediction is applied to the bounding box found by a face detector, as well as to a number of shifted bounding boxes. The bounding box that produces the minimum PLS residual is selected for predicting the pose parameters using the trained regression. This results in a time-consuming estimator since both many bounding boxes must be considered and latent-space PLS projections must be performed at runtime. The advantage of our method with respect to [START_REF] Haj | On partial least squares in head pose estimation: How to simultaneously deal with misalignment[END_REF] is that bounding box shifting is embedded in the learning process such that the optimal shift becomes part of the output, which yields a computationally efficient head pose predictor.

A short version of this paper was recently published [START_REF] Drouard | Head pose estimation via probabilistic high-dimensional regression[END_REF]. The main addition with respect to [START_REF] Drouard | Head pose estimation via probabilistic high-dimensional regression[END_REF] is the use of a partiallylatent response variable associated with the proposed inverse regression method. This can be viewed as an augmented model and its incorporation in our formulation yields excellent results because it enables robustness to phenomena that make the head-pose estimation problem difficult. Both the proposed method (Section III) and the algorithm (Appendix A) are described in more detail than in [START_REF] Drouard | Head pose estimation via probabilistic high-dimensional regression[END_REF], which makes the paper self-contained, self-explanatory, and enables others to easily reproduce the results. Additionally, Section V contains an extended set of results and comparisons, using three publicly available datasets.

III. PARTIALLY-LATENT MIXTURE OF LINEAR REGRESSIONS

In this section we summarize the generative inverse regression method proposed in [START_REF] Deleforge | High-dimensional regression with gaussian mixtures and partially-latent response variables[END_REF]. The method is referred to as Gaussian locally linear mapping (GLLiM). GLLiM interchanges the roles of the input and output variables, such that a low-dimensional to high-dimensional regression problem is solved instead of a high-dimensional to low-dimensional one. The immediate consequence of using such an inverse regression model is a dramatic reduction in the number of model parameters, thus facilitating the task of training.

An additional advantage of this method is that it can be trained by adding a latent part to the output: while the highdimensional input remains fully observed, the low-dimensional output is a concatenation of a multivariate observed variable and a multivariate latent variable. This is referred to as hybrid-GLLiM (or hGLLiM). The latent part of the output variable has a principled mathematical definition (see below) but it does not necessarily have a physical meaning. The main idea (of introducing a partially latent output) relies on the fact that variabilities present in the high-dimensional face descriptor depend on head pose, on face to bounding-box alignment, and on other phenomena (face shapes, facial expressions, gender, hair and skin color, age, etc.) that are not relevant for the task at hand. The latent part of the regression output has the interesting feature of gathering information other than pose parameters and bounding-box parameters.

A. Inverse Regression

Let X, Y be two random variables, such that X ∈ R L denotes the low-dimensional output, e.g. pose parameters and Y ∈ R D (D L) denotes the high-dimensional input, e.g. feature vector. The goal is to predict a response X given both an input Y and the model parameters θ, i.e. p(X|Y ; θ). We consider the inverse low-to-high mapping from the output variable X to the input variable Y . This, possibly nonlinear, mapping is modeled by a mixture of locally-affine transformations:

Y = K k=1 I(Z = k)(A k X + b k + e k ), (1)
where I is the indicator function, and Z is the standard discrete latent variable:

Z = k if and only if Y is the image of X by the affine transformation A k X + b k , with A k ∈ R D×L and b k ∈ R D
, and e k ∈ R D is an error vector capturing both the high-dimensional observation noise and the reconstruction error due to the piecewise approximation. The missing-data variable Z allows one to write the joint probability of X and Y as the following mixture:

p(Y = y, X = x; θ) = K k=1 p(Y = y|X = x, Z = k; θ) × p(X = x|Z = k; θ)p(Z = k; θ), (2)
where θ denotes the model parameters and y and x denote realizations of Y and X respectively. Assuming that e k is a zero-mean Gaussian variable with diagonal covariance matrix Σ k ∈ R D×D with diagonal entries σ k1 , . . . , σ kD , we obtain that p(y|x, Z = k; θ) = N (y;

A k x + b k , Σ k ).
If we further assume that X follows a mixture of Gaussians via the same assignment Z = k, we can write that p(x|Z

= k; θ) = N (x; c k , Γ k ) and p(Z = k; θ) = π k , where c k ∈ R L , Γ k ∈ R L×L and K k=1 π k = 1.
Note that this representation induces a partition of R L into K regions R k , where R k is the region where the transformation (A k , b k ) is most likely invoked, e.g. Fig. 1. This model is described by the parameter set

θ = {c k , Γ k , π k , A k , b k , Σ k } K k=1 . (3)
Notice that the number of model parameters θ is dictated by the number of parameters of a multivariate Gaussian distribution and by the number of Gaussian components (K). However, the number of parameters of an unconstrained GMM is quadratic in the dimension of the variable. Hence, the size of the training dataset required to reliably estimate a conventional GMM would become prohibitively high for the dimension considered in the paper (D = 1888). This is why an inverse regression strategy is adopted, making the number of parameters linear in the input variable dimension rather than quadratic. This drastically reduces the model size in practice, making it tractable.

B. Inverse Regression with Partially Latent Output

We now extend the previous model such that one can train the inverse regression in the presence of partially latent output: hybrid-GLLiM. While the high-dimensional variable Y remains unchanged, i.e. fully observed, the low-dimensional variable is a concatenation of an observed variable T ∈ R Lt and a latent variable W ∈ R Lw , namely X = [T ; W ], where [•; •] denotes vertical vector concatenation and with L t + L w = L. Hybrid-GLLiM can be seen as a latent-variable augmentation of standard regression. It can also be seen as a semi-supervised dimensionality reduction method since the unobserved low-dimensional variable W must be recovered from realizations of the observed variables Y and T .

The decomposition of X implies that some of the model parameters must be decomposed as well, namely c k , Γ k and A k . Assuming the independence of T and W given Z we have

c k = c t k c w k , Γ k = Γ t k 0 0 Γ w k , A k = A t k A w k . (4
) It follows that (1) rewrites as

Y = K k=1 I(Z = k)(A t k T + A w k W + b k + e k ), (5)
While the parameters to be estimated are the same, i.e. (3), there are two sets of missing variables,

Z 1:N = {Z n } N n=1 ∈ {1 . . . K} and W 1:N = {W n } N n=1 ∈ R Lw , associated with the training data (y, t) 1:N = {y n , t n } N
n=1 given the number of K of affine transformations and the latent dimension L w . Also notice that the means {c w k } K k=1 and covariances {Γ w k } K k=1 must be fixed to avoid non-identifiability issues. Indeed, changing their values corresponds to shifting and scaling the latent variables {W n } N n=1 which are compensated by changes in the parameters of the affine transformations {A w k } K k=1 and {b w k } K k=1 . This identifiability problem is the same as the one encountered in latent variable models for dimension reduction and is always solved by fixing these parameters. Following [START_REF] Ghahramani | The EM algorithm for mixtures of factor analyzers[END_REF] and [START_REF] Tipping | Mixtures of probabilistic principal component analyzers[END_REF], the means and covariances are fixed to zero and to the identity matrix respectively: c w k = 0, Γ w k = I, ∀k ∈ {1 . . . K}.

The corresponding EM algorithm consists of estimating the parameter set θ that maximizes

θ (i) = argmax θ E Z [log p((x, y, W , Z) 1:N ; θ| (x, y) 1:N ; θ (i-1) )] . (6)
Using that W 1:N and T 1:N are independent conditionally on Z 1:N and that {c w k } K k=1 and {Γ w k } K k=1 are fixed, maximizing (6) is then equivalent to

θ (i) = argmax θ E r (i) Z [E r (i) W |Z [log p(y 1:N | (t, W, Z) 1:N ; θ)] + log p((t, Z) 1:N ; θ)] , (7) where r (i)
Z and r

(i)
W |Z denote the posterior distributions

r (i) W |Z = p(W 1:N |(y, t, Z) 1:N ; θ (i-1) ), (8) r (i)
Z = p(Z 1:N |(y, t) 1:N ; θ (i-1) ). (9)
It follows that the E-step splits into two steps, an E-W step and an E-Z step. Details of the hybrid-GLLiM algorithm are provided in Appendix A. Note that the model described in Section III-A simply corresponds to L w = 0, hence this algorithm can be used to solve both models (fully observed output and partially observed output).

C. Forward Predictive Distribution

Once the model parameters θ are estimated, one obtains the low-dimensional to high-dimensional inverse predictive distribution as detailed in [START_REF] Deleforge | High-dimensional regression with gaussian mixtures and partially-latent response variables[END_REF]. More interesting, it is also possible to obtain the desired high-dimensional to low-dimensional forward predictive distribution:

p(x|y; θ * ) = K k=1 ν * k N (x; A * k y + b * k , Σ * k ) (10)
with

ν * k = π * k N (y; c * k , Γ * k ) K j=1 π * j N (y; c * j , Γ * j ) (11)
which is also a Gaussian mixture conditioned by the parameters θ * :

θ * = {c * k , Γ * k , π * k , A * k , b * k , Σ * k } K k=1 . (12)
A prominent feature of this model is that the parameters θ * can be expressed analytically from the parameters θ (please consult the analytical expressions of these parameters in Appendix A):

c * k = A k c k + b k , (13)
Γ * k = Σ k + A k Γ k A k , ( 14
) π * k = π k , (15)
A * k = Σ * k A k Σ -1 k , (16) b
* k = Σ * k Γ -1 k c k -A k Σ -1 k b k , (17)
Σ * k = Γ -1 k + A k Σ -1 k A k . -1 (18)
Using [START_REF] Deleforge | High-dimensional regression with gaussian mixtures and partially-latent response variables[END_REF] one can predict an output x corresponding to a test input ŷ x = f (ŷ) with:

f (y) = E [x|y; θ * ] = K k=1 ν * k (A * k y + b * k ) . (19)

IV. IMPLEMENTATION DETAILS

The proposed head-pose estimation method is implemented as follows. We use the Matlab computer vision toolbox implementation of the face detector of [START_REF] Viola | Rapid object detection using a boosted cascade of simple features[END_REF] as we found that this method yields very good face detections and localizations for a wide range of face orientations, including side views. The Matlab implementation of [START_REF] Viola | Rapid object detection using a boosted cascade of simple features[END_REF] offers three different trained classifiers for face detection: two of them for frontal-view detection and one for profile-view detection. These three classifiers yield different results for face detection in terms of bounding-box location and size. The results of face detection using these three classifiers are then combined together for both training and testing of our method. For each face detection the associated bounding box is resized to patches of 64 × 64, converted to a grey-level image to which histogram equalization is then applied. A HOG descriptor is extracted from this resized and histogram-equalized patch. To do so, we build a HOG pyramid (p-HOG) by stacking HOG descriptors at multiple resolutions. The following parameters are used to build p-HOG descriptors: Three HOG descriptors are computed, one for each cell resolution, which are then stacked to form a high-dimensional vector y ∈ R D , with D = 1888.

The dimension of the output variable x ∈ R L depends on the number of pose parameters (up to three angles: yaw, pitch and roll), the bounding-box shift parameters (horizontal and vertical shifts) and the number of latent variables. Hence the output dimension may vary from L = 1 (one angle, no shift, no latent variable) to L = 10 (three angles, two shifts, four latent variables). We used the algorithm detailed in Appendix A to learn the model parameters, either with fully observed output variables (pose angles and bounding-box shifts), or with both observed and latent output variables. From these inverse regression parameters, we derive the forward parameters θ * , i.e. equations ( 13)-( 18) that allow us to estimate the forward predictive distribution [START_REF] Deleforge | High-dimensional regression with gaussian mixtures and partially-latent response variables[END_REF], and hence to predict an output from a test input, i.e. [START_REF] Raytchev | Head pose estimation by nonlinear manifold learning[END_REF]. The joint estimation of the head-pose angles and of the bounding-box shift is achieved iteratively in the following way. The current bounding-box location, u ∈ R 2 , is used to build a feature vector y from which both a head pose x h and a bounding-box shift x b are predicted. The latter is then used to update the bounding-box location, to build an updated feature vector and to predict an updated head pose and a new bounding-box shift. This iterative prediction is described in detail in Algorithm 1.

V. EXPERIMENTS

In this section we present an experimental evaluation of the proposed head-pose estimation methodology. The experiments are carried out with three publicly available datasets: the Prima return head-pose x h and bounding-box location x b 8: end procedure dataset [START_REF] Gourier | Estimating face orientation from robust detection of salient facial features[END_REF], the BIWI dataset [START_REF] Fanelli | Random forests for real time 3d face analysis[END_REF], and the McGill real-world face video dataset [START_REF] Demirkus | Robust semi-automatic head pose labeling for real-world face video sequences[END_REF], [START_REF] Demirkus | Hierarchical temporal graphical model for head pose estimation and subsequent attribute classification in real-world videos[END_REF] ]. Yaw values corresponding to each video frame are estimated using a two-step labeling procedure that provides the most likely angle as well as a degree of confidence. The labeling consists of showing images and possible angle values to human experts, i.e. [START_REF] Demirkus | Robust semi-automatic head pose labeling for real-world face video sequences[END_REF].

We adopted the leave-one-out evaluation protocol at the individual person level. More precisely, all the images/frames associated with one participant are left aside and used for testing, while the remaining data were used for training. As a measure of performance, we evaluated the absolute error between an estimated angle and the ground-truth value, then we computed the mean absolute error (MAE) and standard deviation (STD) over several tests.

We experimented with the following variants of the proposed algorithm (notice that the number of pose angles depends on the dataset, i.e. yaw, pitch and roll (BIWI), yaw and pitch (Prima) or yaw (McGill)):

• GLLiM pose learns and predicts with one, two, or three pose angles; • hGLLiM pose-d learns with the pose parameters as well as with partially latent output, where the number d of latent variables varies between 1 and 4; • GLLiM pose&bb learns and predicts with both pose angles and bounding-box shifts, and

• hGLLiM pose&bb-d learns and predicts with pose angles, bounding-box shifts and partially latent output.

An important aspect of any head-pose method is the way faces are detected in images. We used manually annotated bounding boxes, whenever they are available with the datasets. Otherwise, we used bounding boxes provided with a face detector, e.g. [START_REF] Viola | Rapid object detection using a boosted cascade of simple features[END_REF]. To evaluate the robustness to inaccurate face localization, we introduced random shifts, drawn from a Gaussian distribution, and we used these shifts in conjunction with GLLiM pose&bb and with hGLLiM pose&bb-d to learn the regression parameters and to predict the correct boundingbox location. Notice that in the case of the latter algorithms, the prediction is run iteratively, i.e. Algorithm 1: Extract a HOG vector, predict the pose and the shift, extract a HOG vector from the shifted bounding box, predict the pose and the shift, etc. This is stopped when the shift becomes very small and, as it can be seen below, it considerably improves the quality of the head-pose method. The number K of Gaussian components is an important parameter, which in our model corresponds to the number of affine mappings. We carried out several experiments in order to evaluate the quality of the results obtained by our method as a function of the number of affine transformations in the mixture. We use the GLLiM pose variant of our algorithm with three different face detection options: manual annotation, manual annotation perturbed with additive Gaussian noise, and automatic face detection. We trained these three versions of GLLiM pose with K varying from 1 to 100.

In order to determine the optimal number of affine mappings, K, associated with GLLiM, we use two criteria, the Bayesian information criterion (BIC) which is an informationtheoretic criterion generally used for model selection, and an experimental figure of merit based on the mean absolute error (MAE). We learned several models for different values of K using the Prima dataset. We seek the model that yields low BIC and MAE scores. The BIC and MAE values are plotted as a function of K in Table I and in Fig. 2 and Fig. 3. These curves show the same behavior: as the number of affine mappings increases from K = 1 to approximatively K = 30, both the BIC and MAE scores decrease, then the curve slopes become almost horizontal. Both BIC and MAE reach the lowest score for K = 50. This behavior can be explained as follows. When K < 5 the model is not flexible enough to take into account the apparently non-linear mapping between HOG features and head-pose parameters. It can be observed from Fig. 3 that a large value for K increases the model accuracy. As expected, the computational complexity increases with K as well. Indeed, the number of model parameters is linear in the number of mixture components and hence the size of the training dataset must be increased as well. It is well known that a large number of components in a mixture model presents the risk of overfitting. It is interesting to notice that BIC (derived [START_REF] Fanelli | Random forests for real time 3d face analysis[END_REF] USES DEPTH DATA ONLY AND [START_REF] Wang | Head pose estimation with combined 2D SIFT and 3D HOG features[END_REF] AND [START_REF] Ghiass | Highly accurate and fully automatic head pose estimation from a low quality consumer-level rgb-d sensor[END_REF] from information theory) and MAE (based on experiments with the data) yield the same optimal value, namely K ≈ 50.

The proposed algorithms were compared with the following state-of-the-art head-pose estimation methods: the neuralnetwork based methods of [START_REF] Stiefelhagen | Estimating head pose with neural networks -results on the pointing04 ICPR workshop evaluation data[END_REF], [START_REF] Gourier | Head pose estimation on low resolution images[END_REF] and of [START_REF] Ahn | Real-time head orientation from a monocular camera using deep neural network[END_REF] , the method of [START_REF] Demirkus | Soft biometric trait classification from real-world face videos conditioned on head pose estimation[END_REF] based on dictionary learning, the graphicalmodel method of [START_REF] Demirkus | Probabilistic temporal head pose estimation using a hierarchical graphical model[END_REF], the template based method of [START_REF] Zhu | Face detection, pose estimation, and landmark localization in the wild[END_REF], the supervised non-linear optimization method of [START_REF] Xiong | Supervised descent method and its applications to face alignment[END_REF], the optimization method of [START_REF] Ghiass | Highly accurate and fully automatic head pose estimation from a low quality consumer-level rgb-d sensor[END_REF], and the random-forest methods of [START_REF] Fanelli | Random forests for real time 3d face analysis[END_REF] and of [START_REF] Wang | Head pose estimation with combined 2D SIFT and 3D HOG features[END_REF]. Additionally, we benchmarked the following regression methods: support vector regression (SVR) [START_REF] Smola | A tutorial on support vector regression[END_REF], Gaussian process regression (GPR) [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], and partial least squares (PLS) [START_REF] Abdi | Partial least square regression (PLS regression)[END_REF], as they are widely known and commonly used regression methods for which software packages are publicly available. Notice that some of these methods estimate only one parameter, i.e. the yaw angle [START_REF] Zhu | Face detection, pose estimation, and landmark localization in the wild[END_REF], [START_REF] Demirkus | Soft biometric trait classification from real-world face videos conditioned on head pose estimation[END_REF]- [START_REF] Xiong | Supervised descent method and its applications to face alignment[END_REF], while the random-forest methods of [START_REF] Fanelli | Random forests for real time 3d face analysis[END_REF], [START_REF] Ghiass | Highly accurate and fully automatic head pose estimation from a low quality consumer-level rgb-d sensor[END_REF] and [START_REF] Wang | Head pose estimation with combined 2D SIFT and 3D HOG features[END_REF] use depth information available with the BIWI (Kinect) dataset.

Table II, Table III, and Table IV show the results of headpose estimation obtained the Prima, BIWI, and McGill datasets, respectively. The † symbol indicates that the results are those reported by the authors while the ‡ symbol indicates that the results are obtained using either publicly available software packages or our own implementations. In the case of the Prima dataset, GLLiM pose and hGLLiM pose yield the best results. We note that hGLLiM pose&bb variants of the algorithm (simultaneous prediction of pose, bounding-box shift and partially-latent output) increase the confidence (low STD). Table III shows the results obtained with the BIWI datasets. As already mentioned, [START_REF] Fanelli | Random forests for real time 3d face analysis[END_REF] uses depth information and [START_REF] Wang | Head pose estimation with combined 2D SIFT and 3D HOG features[END_REF], [START_REF] Ghiass | Highly accurate and fully automatic head pose estimation from a low quality consumer-level rgb-d sensor[END_REF] use of depth and color information. Overall, the proposed algorithms compare favorably with [START_REF] Fanelli | Random forests for real time 3d face analysis[END_REF]. hGLLiM pose-4 yields the best MAE for the roll angle, while [START_REF] Ghiass | Highly accurate and fully automatic head pose estimation from a low quality consumer-level rgb-d sensor[END_REF] yields the best MAE for pitch and yaw, but with a high standard deviation. Our algorithm estimates the parameters with the highest confidence (lowest standard deviation). Table IV shows the results obtained with the McGill dataset. The ground-truth yaw values in this dataset are obtained by human experts that must choose among a discrete set of 7 values. Clearly, this is not enough to properly train our algorithms. The method of [START_REF] Demirkus | Probabilistic temporal head pose estimation using a hierarchical graphical model[END_REF] yields the best results in terms of RMSE while hGLLiM pose-2 yields the best results in terms of MAE. Notice that PLS yield the highest confidence in this case.

The results summarized in Table II, Table III and Table IV allow to quantify the variants of the proposed algorithm. With manually annotated bounding boxes, e.g. Tables II and III, there is no notable difference between these variants. Whenever the regression is trained with simulated boundingbox shifts, both GLLiM pose&bb and hGLLiM pose perform better than GLLiM pose. It is interesting to note that in the case of the Prima dataset (Table II), hGLLiM pose-4 (the dimension of the latent part of the output is 4) is the bestperforming variant of GLLiM, while in the case of the BIWI dataset (Table III), hGLLiM Pose&bb-4 is the best performing one. We also experimented with a larger latent dimension without improving the accuracy. Concerning the McGill dataset, neither GLLiM pose&bb nor hGLLiM pose&bb could be trained because the ground-truth face bounding boxes are not available. In the case of this dataset, we used a face detector both for training and for testing. The fact that, overall, hGLLiM pose performs better than GLLiM pose validates the advantage of adding a latent component to the output variable. The latter "absorbs" various phenomena that would otherwise affect the accuracy of the pose parameters.

As already mentioned, GLLiM pose&bb and hGLLiM pose&bb are applied iteratively, until there is no improvement in the predicted output. Fig. 4 plots the MAE of yaw and pitch as a function of the number of iterations of GLLiM pose&bb. As it can be observed from these curves, the MAE decreases after two iterations and then it slightly increases and again it decreases. Therefore, in practice we run two iterations of these algorithms.

Finally, we applied hGLLiM Pose&bb-4 to a set of images that are not contained in any of the three datasets. Fig. 5 shows some of these results. The face detector of [START_REF] Viola | Rapid object detection using a boosted cascade of simple features[END_REF] was used to detect faces (red bounding boxes). The output of hGLLiM Pose&bb-4: the shifted bounding boxes (shown in green) and the estimated pitch, roll, and yaw angles (shown with a green arrow).

VI. CONCLUSIONS

In this paper we proposed a solution to the problem of estimating head pose from the bounding box aligned with a face. Instead of extracting facial landmarks, the method directly maps high-dimensional feature vectors (extracted from faces) onto a low-dimensional manifold. The method relies on learning a mixture of linear regressions. The latter is modeled within the framework of generative methods. More precisely, it is assumed that the high-dimensional feature space is generated from a low-dimensional parameter space. Consequently, an inverse regression strategy is adopted: a lowdimensional to high-dimensional regression is learnt, followed by Bayes inversion.

We experimented with four variants of the proposed algorithm: (i) GLLiM pose, which learns and predicts the pose parameters, (ii) hGLLiM pose, which learns and predicts the pose parameters in the presence of latent variable augmentation of the output, (iii) GLLiM pose&bb, which learns and predicts both the pose parameters and bounding-box shifts, and (iv) hGLLiM pose&bb which combines hGLLiM with GLLiM pose&bb. The experiments and benchmarks, carried out with three publicly available datasets, show that the latentaugmentation variants of the algorithm improve the accuracy of the estimation and perform better than several state-of-theart algorithms.

The methodology presented in this paper has not been tuned for the particular application of head-pose estimation. The algorithms may be applied, with minor modifications, to other high-dimensional to low-dimensional mapping problems, e.g. estimation of human pose from color and depth images. It may also be used as input for gaze estimation or for determining the visual focus of attention, e.g. [START_REF] Massé | Simultaneous estimation of gaze direction and visual focus of attention for multi-person-to-robot interaction[END_REF].

In the future we plan to extend the method to the problem of tracking pose parameters. Indeed, a natural extension of the proposed method is to incorporate a dynamic model to better predict the output variables over time. Such a model could be simultaneously applied to the pose parameters and to the bounding-box shifts. Hence, one can track the image region of interest and the pose parameters in a unified framework.

APPENDIX A EM FOR GLLIM AND HYBRID-GLLIM

This appendix details the EM algorithm that estimates the parameters of the regression method described in Section III. The interest reader is referred to [START_REF] Deleforge | High-dimensional regression with gaussian mixtures and partially-latent response variables[END_REF] for an in-depth description and discussion. Once initialized, at each iteration i, the algorithm alternates between two expectation steps, E-Z and E-W, and two maximization steps, M-GMM and M-mapping:

• E-W-step: Given the current parameter estimates θ (i-1) , the posterior probability r

Σ (i-1) k -1 × y n -A t(i-1) k t n -b (i-1) k . ( 21
)
• E-Z-step: The posterior probability r

for all n and k, where p(y n , t n |Z n = k; θ (i-1) )

= N (t n ; c t k , Γ t k )N (y n ; d k , Φ k ), [START_REF] Sundararajan | Head pose estimation in the wild using approximate view manifolds[END_REF] with:

d k = A t(i-1) k t n + b (i-1) k , Φ k = A w(i-1) k A w(i-1) k + Σ (i-1) k .
The maximization [START_REF] Bay | Speeded-up robust features (surf)[END_REF] can then be performed using the posterior probabilities r correspond to those of a standard Gaussian mixture model on T 1:N , so that we get straightforwardly:

c t(i) k = N n=1 ρ (i) nk t n , (24)
Γ t(i) k = N n=1 ρ (i) nk (t n -c t(i) k )(t n -c t(i) k ) (25)
π (i) k = N n=1 r (i) nk N . ( 26
)
• M-mapping-step: The updating of mapping parameters {A k , b k , Σ k } K k=1 is also in closed-form. The affine transformation matrix is updated with:

A (i) k = Y (i) k X (i) k (S x(i) k + X (i) k X (i) k ) -1 (27)
where:

X (i) k = ρ (i) 1k (x (i) 1k -x (i) k ), . . . , ρ (i)
nk (x

(i) N k -x (i) k ) , Y (i)
k = ρ The intercept parameters are updated with:

b (i) k = N n=1 ρ (i) nk (y n -A (i) k x (i) nk ). (28)
The noise covariance matrices are updated with: where the diag{•} operator sets all the off-diagonal entries to 0. • Initialization: Initial parameters θ (0) are obtained by fitting a GMM with K components to the joint outputinput training dataset {t n , y n } N n=1 .

Σ (i) k = diag A w(i) k S w(i) k A w(i) k + (29)
Note that the following derivations are also valid for the estimation of the parameter set θ in (3), which corresponds to L w = 0, hence the E-W step disappears.

Algorithm 1 : repeat 3 :

 13 Iterative prediction Require: Bounding-box location u and forward model parameters θ * 1: procedure HEADPOSEESTIMATION(u,θ * ) 2Build y from current bounding-box location u 4:Predict x = [x h ; x b ] from y using[START_REF] Raytchev | Head pose estimation by nonlinear manifold learning[END_REF]

5 :

 5 Update the bounding-box location u = u + x b 6: until x b ≤ 7:

Fig. 2 .

 2 Fig. 2. The Bayesian information criterion (BIC) as a function of the number of affine transformations in GLLiM. These experiments use the Prima dataset with the leave-one-out protocol.

Fig. 3 .

 3 Fig. 3. Mean absolute error (MAE) in degrees, for pitch (top) and yaw (bottom), as a function of the number of affine transformations in the mixture of linear regression model. We used GLLiM pose to learn the model parameters independently for pitch and yaw. The three curves correspond to the following face detection cases: manual annotation (red curve), manual annotations with additive noise (blue), and automatic face detection (magenta). These experiments use the Prima dataset with the leave-one-out protocol.

Fig. 4 .

 4 Fig. 4. Mean absolute error (MAE) for pitch and yaw as a function of the number of iterations of GLLiM pose&bb.

Fig. 5 .

 5 Fig. 5. Examples of simultaneous estimation of head-pose angles and of bounding-box shifts. The initial bounding box (found with an automatic face detector) is shown in red. The estimated bounding box is shown in green.

W

  |Z in[START_REF] Osadchy | Synergistic face detection and pose estimation with energy-based models[END_REF] is fully determined by the distributions p(w n |Z n k, t n , y n ; θ (i-1) ) for all n and k, which can be shown to be Gaussian. Their covariance matrices S

  nk = p(Z n = k|t n , y n ; θ (i-1) ) = = π (i-1) k p(y n , t n |Z n = k; θ (i-1)) n , t n |Z n = j; θ (i-1) )

  nk = [t n ; µ w(i) nk ] ∈ R L .The M-step can be divided into two separate steps.• M-GMM-step: The updating of parameters π

Vincent

  Drouard received the M.Sc. in applied mathematics from the school of engineering, university of Nice Sophia-Antipolis in 2013. He is currently a Ph.D. student at the university of Grenoble Alpes (UGA) and member of the PERCEPTION team at INRIA Grenoble Rhône-Alpes. In 2015 he received the best student paper award (second place) at the IEEE International Conference on Image Processing.

  Cell resolutions: 32 × 32, 16 × 16 and 8 × 8, and • Number of orientation bins: 8

• Block resolution: 2 × 2 cells; •

:

  • The Prima head pose dataset consists of 2790 images of 15 persons recorded twice. Pitch values lie in the interval [-60 • , 60 • ], and yaw values lie in the interval [-90 • , 90 • ] with a 15 • step. Thus, there are 93 poses available for each person. All the recordings were achieved with the same background. One interesting feature of this dataset is the the pose space is uniformly sampled.

	The dataset
	is annotated such that a face bounding box (manually
	annotated) and the corresponding yaw and pitch angle
	values are provided for each sample.
	• The Biwi Kinect head pose dataset consists of video
	recordings of 20 people (16 men, 4 women, some of
	them recorded twice) using a Kinect camera. During
	the recordings, the participants freely move their head
	and the corresponding head angles lie in the intervals
	[-60

• , 60 • ] (pitch), [-75 • , 75 • ] (yaw), and [-20 • , 20 • ] (roll). Unlike the Prima dataset, the parameter space is not evenly sampled. The face centers (nose tips) were detected on each frame in the dataset, which allowed us to automatically extract a bounding box for each sample. • The McGill real-world face video dataset consists of 60 videos (a single participant per video, 31 women and 29 men) recorded with the objective of studying unconstrained face classification. The videos were recorded in different environments (both indoor or outdoor) thus resulting in arbitrary illumination conditions and background clutter. Furthermore, the participants were completely free in their behaviors and movements. Yaw angles range in the interval [-90 • , 90 •

TABLE I THE

 I BIC SCORE FOR SEVERAL MODELS LEARNED WITH DIFFERENT VALUES OF K USING THE PRIMA DATASET. WE USE GLLIM POSE TO LEARN EACH MODEL WITH DIFFERENT INPUT DATA (FIG. 2): ANNOTATED FACE POSITION (AFP), ADDING MANUAL NOISE TO THE FACE POSITION (MNA) AND USING A FACE DETECTOR (FD), THE OPTIMAL BIC SCORES ARE IN BOLD.

	Data	K = 1	K = 5	K = 25	K = 50	K = 100
	AFP	-6.0608 -6.5845	-6.822	-6.8429	-6.8173
	MNA	-5.4554 -5.6018 -5.6491 -5.6688	-5.6455
	FD	-4.0596 -4.2602 -4.3144 -4.3366	-4.2307

TABLE II MEAN

 II ABSOLUTE ERROR (MAE) AND STANDARD DEVIATION (STD) (IN DEGREES) OBTAINED WITH VARIOUS HEAD-POSE METHODS, REGRESSION METHODS, AND OUR METHOD USING THE PRIMA DATASET. THIS DATASET CONTAINS MANUALLY ANNOTATED BOUNDING BOXES OF FACES AND THE CORRESPONDING PITCH AND YAW ANGLES. IN ORDER TO TEST THE ROBUSTNESS WE SIMULATED SHIFTED BOUNDING BOXES. THE BEST RESULTS ARE IN BOLD.

		Manually annotated bounding boxes	Manual annotation + simulated shifts
		Pitch		Yaw	Pitch	Yaw
	Method	MAE	STD	MAE	STD	MAE	STD	MAE	STD
	Stiefelhagen [37] ‡	9.7	-	9.5	-	-	-	-	-
	Gourier et al. [38] ‡	12.1	-	7.3	-	-	-	-	-
	GPR [39] †	11.94	10.19	15.04	12.24	19.96	16.58	23.69	18.16
	PLS [40] †	12.25	9.73	13.38	10.8	17.77	14.47	17.34	13.94
	SVR [41] †	11.25	9.42	12.82	10.99	17.09	14.81	17.27	14.09
	GLLiM pose	8.41	10.65	7.87	8.08	15.99	16.69	13.66	14.78
	hGLLiM pose-2	8.47	10.35	7.93	7.9	12.64	14.49	11.51	11.37
	hGLLiM pose-4	8.5	10.8	7.85	7.98	12.03	13.92	10.78	9.77
	GLLiM pose&bb	-	-	-	-	13.13	13.65	11.3	10.55
	hGLLiM pose&bb-2	-	-	-	-	12.52	12.44	11.04	9.7
	hGLLiM pose&bb-4	-	-	-	-	12.12	12.85	11.27	9.53

TABLE III THE

 III MEAN ABSOLUTE ERROR (MAE) AND STANDARD DEVIATION (STD), EXPRESSED IN DEGREES, OBTAINED WITH VARIOUS HEAD-POSE METHODS, REGRESSION METHODS, AND OUR METHOD USING THE BIWI DATASET. THIS DATASET CONTAINS ANNOTATED BOUNDING BOXES OF FACES AND THE CORRESPONDING PITCH, YAW, AND ROLL ANGLE VALUES. IN ORDER TO TEST THE ROBUSTNESS WE SIMULATED SHIFTED BOUNDING BOXES BOTH FOR TRAINING AND FOR TESTING. THE BEST RESULTS ARE IN BOLD. NOTE THAT

  USE BOTH COLOR AND DEPTH INFORMATION. PAPERS USING DEPTH DATA ARE MARKED WITH A * .

			Manually annotated bounding boxes			Manual annotation + simulated shifts
		Pitch	Yaw	Roll	Pitch	Yaw	Roll
	Method	MAE	STD	MAE	STD	MAE	STD	MAE	STD	MAE	STD	MAE	STD
	Ghiass et al. [17] ‡ *	0.1	6.7	0.2	8.7	0.3	9.3	-	-	-	-	-	-
	Fanelli et al. [15] ‡ *	3.8	6.5	3.5	6.8	5.4	6.0	-	-	-	-	-	-
	Wang et al. [42] ‡ *	8.5	11.1	8.8	14.3	7.4	10.8	-	-	-	-	-	-
	GPR [39] †	9.64	8.85	7.72	7.17	6.01	6.29	10.77	9.45	9.06	8.33	6.54	6.72
	PLS [40] †	7.87	6.73	7.35	6.06	6.11	5.9	10.32	8.64	8.67	7.7	6.69	6.74
	SVR [41] †	7.77	6.85	6.98	6.26	5.14	5.96	10.82	9.22	9.14	8.32	6.26	7.16
	GLLiM pose	5.77	5.77	4.48	4.33	4.71	5.31	11.33 11.58	10.2	11.34	7.76	8.02
	hGLLiM pose-2	5.57	5.48	4.33	4.68	4.37	5.09	9.04	9.13	7.65	8.3	6.3	6.75
	hGLLiM pose-4	5.43	5.44	4.24	5.37	4.13	4.86	8.45	8.41	6.93	7.72	6.12	6.8
	GLLiM pose&bb	-	-	-	-	-	-	8.49	8.79	6.86	7.3	6.57	6.95
	hGLLiM pose&bb-2	-	-	-	-	-	-	7.81	7.68	6.41	7.19	5.75	6.68
	hGLLiM pose&bb-4	-	-	-	-	-	-	7.65	8.0	6.06	6.91	5.62	6.35

TABLE IV ROOT

 IV MEAN SQUARE ERROR (RMSE), MEAN ABSOLUTE ERROR (MAE) AND STANDARD DEVIATION (STD) (IN DEGREES) OBTAINED WITH VARIOUS HEAD-POSE METHODS, REGRESSION METHODS, AND OUR METHOD USING THE MCGILL REAL-WORLD DATASET. THIS DATASET CONTAINS ANNOTATED YAW ANGLES. BOUNDING BOXES ARE LOCATED WITH A FACE DETECTOR. THE BEST RESULTS ARE IN BOLD.

		Bounding boxes based on face detection
			Yaw
	Method	RMSE	MAE	STD
	Demirkus et al. [43] ‡	> 40	-	-
	Xiong and De la Torre [45] ‡	29.81	-	-
	Zhu and Ramanan [28] ‡	35.70	-	-
	Demirkus et al. [44] ‡	12.41	-	-
	GPR [39] †	23.18	16.22	16.71
	PLS [40] †	22.46	15.56	16.2
	SVR [41] †	21.13	15.25	18.43
	GLLiM pose	26.62	13.1	23.17
	hGLLiM pose-2	24.0	11.99	20.79
	hGLLiM pose-4	24.25	12.01	21.06

Supplementary material, that includes a complete Matlab package and examples of trained regression functions, is publicly available at https://team. inria.fr/perception/research/head-pose.


Joanna Daaboul
email: joanna.daaboul@utc.fr

Julien Le Duigou
email: julien.le-duigou@utc.fr

Diana Penciuc
email: diana.penciuc@utc.fr

Benoît Eynard
email: benoit.eynard@utc.fr

Reverse Logistics: network design based on life cycle assessment

Keywords: Reversed logistics design, LCA, wrought aluminum

When aiming to a more sustainable world, enterprises such as aircraft and automobile industries are highly interested in light weight components and solutions. Of these solutions are aluminum wrought alloys that offer high potentials for dramatic weight reduction of structural parts. Nevertheless, the production of virgin aluminum is, however, highly energy consuming. Hence, and in SuPLight project, we are interested in recycled aluminum. The aim of this project is to address new industrial models for sustainable light weight solutionswith recycling in high-end structural components based on wrought alloys. In this article we address the issue of designing the reverse logistics chain assuring the needed volume of recycled aluminum for the production of Lshaped Front Lower Control Arms for personal cars.

Introduction

Until recently, Reverse Logistics (RL) was not given a great deal of attention in organizations. Actually, implementing RL programs to reduce, reuse, and recycle wastes from distribution and other processes generates tangible and intangible value and can lead to better corporate image [START_REF] Carter | Reverse logistics: a review of the literature and framework for investigation[END_REF]. Its main drivers are legislations and directives, consumer awareness and social responsibilities towards environment [START_REF] Ferrer | The impact of remanufacturing in the economy[END_REF][START_REF] Castell | Extended producer responsibility policy in the European Union-a horse or a camel[END_REF][START_REF] Ravi | Analysis of interactions among the barriers of reverse logistics[END_REF].

Another motivating driver of RL is economic factors. RL can generate profits by reselling valuable components or products [START_REF] Toffel | The growing strategic importance of end-of-life product management[END_REF]. Reverse logistics is of high importance for aluminum based products, since the production of virgin aluminum is highly energy consuming. RL operations and chains they support are significantly more complex than traditional manufacturing supply chains [START_REF] Amini | Designing a reverse logistics operation for short cycle time repair services[END_REF][START_REF] Rogers | Going backwards: reverse logistics trends and practices[END_REF][START_REF] Mitra | Revenue management for remanufactured products[END_REF]. Therefore, and similarly to how companies develop efficient logistics processes for new goods, it is necessary to plan operations for returned goods, taking into consideration that the processes are most probably quite different from those defined for forward distribution [START_REF] Stock | Many happy (product) returns[END_REF]. RL is not a symmetric picture of forward distribution [START_REF] Fleischmann | Quantitative models for reverse logistics: a review[END_REF] and requires different management and planning approaches. For example, it is difficult, in the case of RL, to estimate supply-related parameters such as the unit operational costs directly from reported statistical data.

In this article we are interested in designing an effective reverse logistics network for the L-shaped Front Lower Control Arms (FLCA) (Figure 1) for personal cars [START_REF]Document Deliverable 1.5: SuPLight application scenarios[END_REF] within the SuPLight1 project which is a multidisciplinary research project, combining physics at the atomic scale level, metallurgy, continuum mechanics, structural mechanics, optimization algorithms, tolerance analysis, life cycle analysis, manufacturing and business modeling. The project addresses new industrial models for sustainable light weight solutionswith 75% recycling in high-end structural components based on aluminum wrought alloys.

Figure 1: Front Lower Control Arm (FLCA) [START_REF] Amini | Designing a reverse logistics operation for short cycle time repair services[END_REF] In our case study we are interested in designing the reverse logistics with the following characteristics: 1. The RL network should be environmental friendly as much as possible 2. The different facilities (collection locations, remanufacturing facilities, etc) are not necessarily owned by one company 3. The reverse flow has different sources and might have depending on these sources different disposal routes (reselling, remanufacturing, recycling, etc).

Related Works

Designing a reverse logistics network is usually achieved via mathematical programming. Generally a mixed-integer linear programming model is generated in order to define the optimum collection locations and recycling factories [START_REF] Ramezani | A new multi-objective stochastic model for a forward/reverse logistic network design with responsiveness and quality level[END_REF]. After analyzing the literature we concluded the following:

The literature presents three main types of networks: Forward logistics, reverse logistics, and forward/reverse logistics [START_REF] Ramezani | A new multi-objective stochastic model for a forward/reverse logistic network design with responsiveness and quality level[END_REF]. We are only considering the reverse logistics network.

The most considered decision criterion is the total cost of the network [START_REF] Min | The spatial and temporal consolidation of returned products in a closed-loop supply chain network[END_REF][START_REF] Üster | Benders decomposition with alternative multiple cuts for a multi-product closed-loop supply chain network design model[END_REF][START_REF] Listes | A stochastic approach to a case study for product recovery network design[END_REF][START_REF] Demirel | A mixed-integer programming model for remanufacturing in reverse logistics environment[END_REF][START_REF] Pishvaee | Reverse logistics network design using simulated annealing[END_REF], followed by the service level [START_REF] Du | A bi-objective reverse logistics network analysis for post-sale service[END_REF], and the generated profit [START_REF] Ramezani | A new multi-objective stochastic model for a forward/reverse logistic network design with responsiveness and quality level[END_REF]. In the recent literature review in [START_REF] Ramezani | A new multi-objective stochastic model for a forward/reverse logistic network design with responsiveness and quality level[END_REF], the only objectives considered in RL or forward/reverse chain design are: cost, profit, responsiveness (Service level), source balance, and quality. [START_REF] Zhou | Generic Model of Reverse Logistics Network Design[END_REF] propose a generic model for reverse logistics design considering only costs as decision criteria. [START_REF] Lee | Dynamic network design for reverse logistics operations under uncertainty[END_REF] propose a two-stage stochastic programming model for multiperiod reverse logistics network design with as main decision criteria the different investment and operational costs. [START_REF] Lambert | A reverse logistics decisions conceptual framework[END_REF] consider the environmental regulations as a constraint to their decision model but not as decision criteria. [START_REF] Cardoso | Design and planning of supply chains with integration of reverse logistics activities under demand uncertainty[END_REF] also proposed a mixed integer linear programming model for supply chain planning including reverse logistics activities. Yet they consider the expected net present value (ENPV), as the decision criterion. In our case, we are interested not only in the total cost of the network, but also in the environmental impact of the network as a main decision criterion. Therefore, the remanufacturing and recycling processes impacts on environmental performance of the network need to be considered.

The main decisions of the models are: location/allocation of facilities and transportation values [START_REF] Min | The spatial and temporal consolidation of returned products in a closed-loop supply chain network[END_REF][START_REF] Üster | Benders decomposition with alternative multiple cuts for a multi-product closed-loop supply chain network design model[END_REF][START_REF] Listes | A stochastic approach to a case study for product recovery network design[END_REF][START_REF] Demirel | A mixed-integer programming model for remanufacturing in reverse logistics environment[END_REF][START_REF] Pishvaee | Reverse logistics network design using simulated annealing[END_REF][START_REF] Du | A bi-objective reverse logistics network analysis for post-sale service[END_REF]. Moreover, most of the works found in the literature consider the network design problem with collection, sorting, and disposal facilities owned by one company which is collecting the reverse flow and reusing it. But in many cases, reverse logistics do not only include facilities owned by the company itself, but facilities owned by other partners in the chain. In this case, designing the best RL network consists of 2 main problems: choosing the partner and optimizing the whole network formed of facilities owned by different partners. Since reverse logistics are a part of green supply chains and since we are interested in the environmental performance of the network, the choice of the partners is to be made based on their impact not only on the total cost of the network, but also on the environmental performance of the network.

Most of the works found in the literature consider mainly the remanufacturing option of the disposal step [START_REF] Min | The spatial and temporal consolidation of returned products in a closed-loop supply chain network[END_REF], [START_REF] Üster | Benders decomposition with alternative multiple cuts for a multi-product closed-loop supply chain network design model[END_REF], [START_REF] Demirel | A mixed-integer programming model for remanufacturing in reverse logistics environment[END_REF][START_REF] Pishvaee | Reverse logistics network design using simulated annealing[END_REF][START_REF] Du | A bi-objective reverse logistics network analysis for post-sale service[END_REF] and only few consider the recycling option of the disposal step such as [START_REF] Listes | A stochastic approach to a case study for product recovery network design[END_REF], [START_REF] Wang | A closed-loop logistic model with a spanning-tree based genetic algorithm[END_REF]. But in our case, we are interested in recycling the reverse flow as well as remanufacturing it. In other words, we are interested in different sources of reverse flows which have each different disposal route.

3

A Framework for the design of reverse logistics

Network design method

Considering our analysis of the literature and based on our specific needs, a life cycle approach seems adequate to analyze the performance of the chosen network. It permits assessing the environmental performance of the RL network alone and also its impact on the whole life cycle assessment of the product. Thus, we propose a comparison based method formed of 5 main steps described in Figure 2. Step 1: Context definition In this step, the company and the product (s) are described.

Step 2: Parameters definition Since we are interested in defining the different partners in the RL network as well as the environmental performance of the different processes in the network, a framework linking the processes to the partners was developed.

According to the literature, reverse logistics is formed of 4 main phases: gatekeeping, collection, sorting, and disposal [START_REF] Wang | A closed-loop logistic model with a spanning-tree based genetic algorithm[END_REF][START_REF] Rogers | Going backwards: Reverse logistics trends and practices. Reverse logistics executive council[END_REF][START_REF] Schwartz | Reverse logistics strengthens supply chains[END_REF][START_REF] Stock | Product returns/reverse logistics on Warehousing[END_REF]. The disposal step includes landfilling, remanufacturing, recycling and reselling. For each phase the set of activities to be achieved was identified based on literature analysis. Nevertheless, these activities differ between the design stage and the implementation (execution) stage. Thus, the differentiation between these activities was made and a framework based on design and implementation stages is proposed and presented in Figure 3. The parameters definition consists of identifying the different roles in an RL chain, the responsibilities and activities per role, the requirements for these activities, and the performance indicators per requirement per activity. Four roles were identified: collector, sorter, disposer and the end collector of the reversed flow. A same partner may have more than one role. For each of these roles, the responsibilities and activities were identified for design and implementation stages based on the previously presented framework. These activities are at a macro level and need to be detailed at a micro level but for specific cases or products.

The requirements per activity are case dependent. Most of it depends on the type of the product. The recycling processes for example differ based on the product type. Whereas the general performance indicators are: Total cost of the network including cost of transportation, facilities, purchasing, activities, and stock, Total emission of CO2, an aggregation of main used life cycle assessment indicators such as carbon footprint of product, energy related indicators, chemicals used, PFC/GHG emissions, fresh water consumption, use of renewable/nonrenewable energy, and recycling related indicators such as the amount of recycled material (in our case it is aluminum), the collection quota CQ (the quantity of secondary material, which is recovered by collection systems, related to the total quantity of used products), and the technical recycling quota RQ (it is the relation between the remelted and the collected quantity and describes the yield of technical processes) [START_REF] Quinkertz | A scenario to optimise the energy demand of aluminium production depending on the recycling quota[END_REF]. Moreover, in step 2, case dependent parameters are identified. These are the candidate companies for playing one or more roles in the RL network. For each company a description sheet is filled contating information on the company, some performance indicators, and constraints such as minimum quantity of reverse flow, replenishement lead time, etc.

Step 3: Scenarios definition

The scenario definition is achieved by the expert user in one of 3 ways:

1. Defining any RL network scenario, by choosing for each role one or more companies as well as the quantity of reverse flow and then by choosing the routes between these companies. 2. Defining the RL network step by step, using decision criteria per role. In this case, and for every role, the user may ask to identify the best company based on many decision criteria, such as distance, cost, processes used, labels, reverse flow quality, etc. This option is based on the Analytical Hierarchy Process (AHP). The steps that need to be followed in this option are described in figure 4. AHP is the most used decision support system for supplier selection [START_REF] Van Der Rhee | Understanding trade-offs in the supplier selection process: The role of flexibility, delivery, and value-added services/support[END_REF]. This technique which was developed by Thomas L.

Saaty [START_REF] Saaty | The Analytic Hierarchy Process[END_REF] relies on the expertise of the user in order to generate the weights of attributes. It is based on the comparison of pairs of options and criteria. It has found widespread application in decision-making problems, involving multiple criteria in systems of many levels [START_REF] Liu | The voting analytic hierarchy process method for selecting supplier[END_REF]. Its main advantages are: o consideration of non-tangible subjective attributes, o ability to structure a complex problem, multi-criteria, multi-person and multi-period hierarchically, o capability to investigate each level of the hierarchy separately, and to combine the results as the analysis progresses. o binary comparison of elements (alternatives, criteria and subcriteria), o ease of its IT support. Defining the optimal network based on the most used objective: the total cost using the optimisation option of the prototype. This option is a work-inprogress. The user may choose many scenarios as desired and compare their performance in steps 4 & 5. In other words, the user may compare on the basis of all environmental and economic criterias the optimal scenario (based on cost only) with scenarios built with the decision aid option and with any other scenario built by the user alone without any aid of the system.

Steps 4 & 5: assessment of RL scenarios and results analysis

After one or more scenarios are chosen, the different performance indicators are calculated. A comparision between the scenario is provided. The evaluation of the RL network is achieved via a life cycle assessment software (SimaPro) [START_REF] Suplight | Deliverable 4.3: LCA/LCC tools for lightweight solutions[END_REF]. This is possible within the Suplight project which aims at providing collaborative platform [START_REF] Suplight | Deliverable 3.2: Generic framework for simulation-based optimization[END_REF] allowing the communication and exchange of information between different software of which the Life Cycle Assessment software (LCA) and the RL demonstrator. After the life cycle assesment of the different RL scenarios is received from the LCA software, a full analysis combing this assessment and the one achieved via the RL prototype (total cost of the network, quantity of recycled material, etc.) is provided to the user. A comparison of all created RL scenarios is then delivered to the user. The proposed framework was the basis of the development of the software prototype illustrating our method. It has been designed with a flexible architecture in order to allow portability and interoperability with any external application. The software prototype implementation is based on the MVC software architecture pattern, using HTML5 and CSS to create the view and for the application logic the JAX-RS Java API for RESTful web services which are controlled by JavaScript (JQuery).

The RL prototype is focused on RL scenario definition allowing the user to: choose the partners of the RL chain and their roles and visualize them on a map; define RL routes -a single route is characterized by flow and transportation data (a flow type and quantity exchanged between two partners, transportation mean and distance)-; visualize the summary of the created scenario and display LCA indicators for the current scenario. The total cost is computed for each route. Google services are used for automatic distance calculation and map display based on data provided by the user. The prototype main page and the route definition interface are presented in Figures 5 and6 respectively. On the main page the prototype shows the RL scenario on a map with a description of the different routes including the nodes names (companies' names) with the type of reverse flow, quantity of flow, transportation mode and distance between the two nodes. The different enterprises or facilities chosen for the RL chain are presented on the map. This map is dynamic and shows the chosen facility locations or the chosen suppliers as well as all related costs. This map serves only as a visualization tool for the chosen scenario. This prototype is not intended for advanced planning of a supply chain such as APS software (Advanced Planning and Scheduling). But its main function is to provide a support for decision making in designing a reverse logistics chain. Even though APS software may provide "What if" analysis, it doesn't offer a decision aid method such as AHP, nor does it provide the optimal RL chain among all possible scenarios; it permits only to compare between manually created scenarios. The use of the proposed demonstrator doesn't eliminate the need to use APS software which is necessary to manage a RL network. The demonstrator may be considered as a plug-in for APS software.

Conclusion and perspectives

In this article we presented a framework for designing an RL network based on life cycle assessment. We are interested in identifying the partners of the network and the routes connecting them based on the total cost of the network and its environmental performance as well. In the perspectives of the presented work is the achievement of the optimization option of the prototype as well as applying the method and testing the prototype in the LFCA case study in the Suplight project.

Figure 2 :

 2 Figure 2: Proposed method for reverse logistics network design

  and the composition of the reverse flow Identifying returns within the acceptable defined characteristics and composition Sorting Defining for every return the best fitting sorting option Disposal Sending the returns to their disposal destination Disposal of returns (remanufacturing/ reselling/ recyling/ etc…) Analyzing the returns and define possible sorting routes Sending the collected volume of returns to the sorting

Figure 3 :

 3 Figure 3: Reverse Logistics Framework

Figure 4 :

 4 Figure 4: Steps to use the multi-criteria decision aid option 3.Defining the optimal network based on the most used objective: the total cost using the optimisation option of the prototype. This option is a work-inprogress. The user may choose many scenarios as desired and compare their performance in steps 4 & 5. In other words, the user may compare on the basis of all environmental and economic criterias the optimal scenario (based on cost only) with scenarios built with the decision aid option and with any other scenario built by the user alone without any aid of the system.

Figure 5 :

 5 Figure 5: RL Prototype main page

Figure 6 :

 6 Figure 6: Route definition interface of the RL prototype

http://www.suplight-eu.org/

Acknowledgments

The authors would like to acknowledge the European Commission and the Picardie French Region for their financial support through the SuPlight FP7 Project (grant agreement n° 263302). We also wish to express our gratitude and appreciation to all the project partners for their contribution during the development of various ideas and concepts presented in this paper.