Mostafa Sadeghi

Sylvain Guy

Adrien Raison

Xavier Alameda-Pineda

Radu Horaud
email: radu.horaud@inria.fr

Unsupervised performance analysis of 3D face alignment

We address the problem of analyzing the performance of 3D face alignment (3DFA) algorithms. Traditionally, performance analysis relies on carefully annotated datasets. Here, these annotations correspond to the 3D coordinates of a set of pre-defined facial landmarks. However, this annotation process, be it manual or automatic, is rarely error-free, which strongly biases the analysis. In contrast, we propose a fully unsupervised methodology based on robust statistics and a parametric confidence test. We revisit the problem of robust estimation of the rigid transformation between two point sets and we describe two algorithms, one based on a mixture between a Gaussian and a uniform distribution, and another one based on the generalized Student's t-distribution. We show that these methods are robust to up to 50% outliers, which makes them suitable for mapping a face, from an unknown pose to a frontal pose, in the presence of facial expressions and occlusions. Using these methods in conjunction with large datasets of face images, we build a statistical frontal facial model and an associated parametric confidence metric, eventually used for performance analysis. We empirically show that the proposed pipeline is neither method-biased nor data-biased, and that it can be used to assess both the performance of 3DFA algorithms and the accuracy of annotations of face datasets.

Introduction

The problem of face alignment (FA) is the problem of facial landmark detection and localization from a single RGB image. Face alignment is an important research topic as it provides input to a variety of computer vision tasks, such as head-pose estimation and tracking, face recognition, facial expression understanding, visual speech recognition, etc., [START_REF] Escalera | Special issue on the computational face[END_REF][START_REF] Loy | Special issue on deep learning for face analysis[END_REF]. 2D face alignment (2DFA) has been extensively studied for the last decades, yielding a plethora of methods and algorithms [START_REF] Wu | Facial landmark detection: A literature survey[END_REF]. State of the art 2DFA based on deep neural networks (DNNs) are the best-performing methods in terms of accuracy, invariance with respect to facial appearances, shapes, expressions, as well as in terms of repeatability and reproducibility in the presence of image noise, image resolution, motion blur, lighting conditions and varying backgrounds.

Nevertheless, 2DFA methods yield poor landmark detection and localization performance in the presence of occlusions which occur in case of large poses induced by out-of-imageplane head rotations (self occlusions) as well as by the presence of various objects in the camera field of view, such as glasses, hair, hands and handheld objects, etc. Robust facial landmark detection and localization in the presence of occlusions can only be achieved on the premise that 3D information is taken into account. It is well established that 2D facial landmarks (and, more generally, face images) embed 3D information. This information can be retrieved by fitting a 3D face model to a 2D face image, even if the latter is only partially visible. The process of fitting a 3D model to a 2D image constitutes the basis of training 3D face alignment (3DFA) algorithms.

Consider for example a 3D face model that is parameterized both by identities and by facial deformations, e.g. the parametric 3D deformable model (3DMM) [START_REF] Blanz | A morphable model for the synthesis of 3D faces[END_REF]. The task of fitting 3DMM to an RGB image of a face consists of estimating the parameters of the mapping from the 3D generic model to a particular face, namely the identity and expression parameters, as well as the pose parameters (scale, rotation, translation and projection), e.g. [START_REF] Gou | Shape augmented regression for 3D face alignment[END_REF][START_REF] Zhu | Face alignment across large poses: a 3D solution[END_REF]. Once an optimal set of parameters is found, one can associate 3DMM vertices with facial landmarks. This stays (a) Analysis of 3D landmarks extracted with [START_REF] Bulat | Two-stage convolutional part heatmap regression for the 1st 3D face alignment in the wild (3DFAW) challenge[END_REF] (b) Analysis of ground-truth 3D landmarks from the AFLW2000-3D dataset [START_REF] Zhu | Face alignment across large poses: a 3D solution[END_REF] Fig. 1: Two examples from the AFLW2000-3D dataset (Zhu et al, 2016) (left). The landmarks are mapped onto a statistical frontal landmark model (right) built using the YawDD dataset [START_REF] Abtahi | Yawdd: A yawning detection dataset[END_REF] and [START_REF] Feng | Joint 3D face reconstruction and dense alignment with position map regression network[END_REF], which enables us to verify whether the mapped landmarks lie within their associated ellipsoidal confidence volumes or not.

at the basis of many automatic and semi-automatic methods for annotating 2D faces with 3D landmarks, e.g. [START_REF] Deng | The Menpo benchmark for multi-pose 2D and 3D facial landmark localisation and tracking[END_REF].

Nevertheless, the fitting task just mentioned is a difficult nonlinear optimization problem, in particular in the presence of large poses and of occlusions. In the recent past, a number of methods has been developed to perform this 3D-to-2D fitting process necessary for 3D facial landmark annotation. The performance of the vast majority of existing 3DFA methods inherently rely on the quality of landmark annotation. This is true for training using modern discriminative deep learning methods, but it is true for testing as well. Indeed, to date, algorithm performance is computed empirically by measuring the error between the predicted output and the corresponding groundtruth, e.g. [START_REF] Jeni | The first 3D face alignment in the wild (3DFAW) challenge[END_REF][START_REF] Deng | The Menpo benchmark for multi-pose 2D and 3D facial landmark localisation and tracking[END_REF]. Under these circumstances, annotation errors are likely to bias both parameter estimation (training) and performance evaluation (testing).

There is a lack of a benchmarking methodology that could assess quantitatively and in a completely unsupervised manner the robustness and effectiveness of 3DFA algorithms, namely a method that computes a confidence score that measures algorithm performance in the absence of the ground truth. This is also crucial in order to decide without human intervention, whether a 3DFA method, that is applied to an unknown image of a face with no annotation available, yields an output that is accurate enough to be further used by other algorithms, such as head-gaze estimation, facial expression analysis or lip reading. This paper proposes a methodological framework for assessing the performance of 3DFA algorithms based on robust probability distribution functions and on a statistical confidence test.

Unlike supervised metrics currently in use for 3DFA performance evaluation and based on annotated datasets, the proposed method is fully unsupervised. We show that the robust estimation of the rigid mapping between two sets of 3D facial landmarks, namely (i) a predicted set, associated with a face with unknown identity, pose and expression, and (ii) a model set associated with a statistical frontal face, provides a reliable way to separate face pose (due to head motions) from non-rigid face deformations (due to facial expressions), all in the presence of badly located landmarks.

Using a 3DFA algorithm and a very large and unannotated dataset of face images with large variabilities in pose, expression and identity, we make use of the robust rigid-mapping methodology to build a statistical frontal landmark model and a parameterized confidence score. Based on this pipeline, the proposed performance evaluation protocol proceeds as follows. First, 3D landmarks are extracted from a face image using a 3DFA method. Second, the landmarks are rigidly mapped onto the frontal model. Third, a confidence score is computed for each mapped landmark, thus allowing to assess whether the landmark lies within a confidence volume or not.

We describe in detail an experimental evaluation framework that uses publicly available datasets and 3DFA software packages associated with three published articles and one unpublished paper. We empirically show that our methodological pipeline is neither dataset-nor method-biased. We also show that the proposed framework can be used not only to assess quantitatively the performance of 3DFA algorithms, but also to test the accuracy of automatic and semi-automatic methods currently used for the annotation of face datasets.

The methodology proposed in this paper is illustrated in Fig. 1. The two images (left) are from the AFLW2000-3D dataset [START_REF] Zhu | Face alignment across large poses: a 3D solution[END_REF]. The statistical frontal landmark model (right) is built using the 3DFA method of [START_REF] Feng | Joint 3D face reconstruction and dense alignment with position map regression network[END_REF] and the YawDD dataset [START_REF] Abtahi | Yawdd: A yawning detection dataset[END_REF]. This model characterizes each landmark with an ellipsoidal confidence volume centered at a posterior mean. Fig. 1(a): Landmarks extracted using (Bulat and Tzimiropoulos, 2016) (left) and mapped onto the statistical model (right). In this case, most of the landmarks lie inside their confidence volumes, thus assessing their correctness. Fig. 1(b): Ground-truth landmarks obtained with a semiautomatic annotation process [START_REF] Zhu | Face alignment across large poses: a 3D solution[END_REF] and mapped onto the statistical model (right). One may notice that in this case, many mapped landmarks fall outside their confidence volumes. The benefit of the proposed method is twofold: (i) an unsupervised assessment of the quality of the detected landmarks, and (ii) a robust and expression-and identity-preserving landmark mapping from an arbitrary pose to a frontal pose.

The remainder of this paper is organized as follows. Section 2 reviews the related work. Section 3 summarizes the problem formulation. Section 4 briefly reviews maximum likelihood estimation (MLE) for rigid mapping and describes two robust methods. Section 5 empirically analyses the proposed rigid-mapping methods. Section 6 proposes a methodological pipeline for building a statistical face model and an associated parametric confidence metric. Section 7 presents extensive ex-perimental results, and Section 8 draws some conclusions. 1

Related Work

It is interesting to note that the recently proposed methods for 3DFA lie at the crossroads of deformable shape models, model-based image analysis and neural networks. In order to discuss these links, we introduce some mathematical notations and concepts. Let vector p ∈ P ⊂ R K denote the ensemble of parameters of a 3D face model S (identity, expression and pose), where P is the parameter vector space, K is the number of parameters, and let b ∈ B ⊂ R I×J denote the image of a face from a set of images of size I × J. One class of 3DFA methods directly learns a mapping p = f (b) from a training dataset of M face images and their associated model parameters {b m , p m } M m=1 , e.g. [START_REF] Zhu | Face alignment across large poses: a 3D solution[END_REF][START_REF] Jourabloo | Pose-invariant face alignment via CNN-based dense 3D model fitting[END_REF][START_REF] Feng | Joint 3D face reconstruction and dense alignment with position map regression network[END_REF][START_REF] Deng | The Menpo benchmark for multi-pose 2D and 3D facial landmark localisation and tracking[END_REF].

Another class of methods proceeds in two steps. First, 2D landmarks are extracted from a face image by learning an image-to-landmark mapping u = g(b), from a face image to a set of 2D landmarks u = {v n } N n=1 ∈ N 2×N , and using a training

dataset {b m , u m } M m=1 . Second, a 2D-to-3D mapping s = h(u) is estimated, where s = {x n } N n=1 ∈ R 3×N is a set of N 3D land- marks.
This mapping can be obtained either by learning, using a training dataset conditioned by a 3D model parameterized by p, i.e. {u m , s( p m )} M m=1 , e.g. [START_REF] Zhao | Fast and precise face alignment and 3D shape reconstruction from a single 2D image[END_REF]Bulat andTzimiropoulos, 2016, 2017), or by direct optimization over p of a function that maps a 3D model onto the 2D landmarks, e.g. [START_REF] Yu | Learning dense facial correspondences in unconstrained images[END_REF]. These 3DFA DNN-based methods use a variety of architectures in order to learn the regression functions f (•), g(•) and h(•) mentioned above. Given this variety, it is difficult to directly compare them and assess their merits based on the analysis of the underlying DNN concepts and methodologies. Hence, 3DFA algorithm performance should be measured empirically, as is often the case in deep learning.

To date and to the the best of our knowledge, there has been a handful of 3DFA benchmark datasets and corresponding evaluation metrics, [START_REF] Jeni | The first 3D face alignment in the wild (3DFAW) challenge[END_REF][START_REF] Deng | The Menpo benchmark for multi-pose 2D and 3D facial landmark localisation and tracking[END_REF][START_REF] Sanyal | Learning to regress 3D face shape and expression from an image without 3D supervision[END_REF]. In [START_REF] Jeni | The first 3D face alignment in the wild (3DFAW) challenge[END_REF], four datasets were specifically gathered, annotated and prepared, and two performance metrics were used for this challenge. The BU-4DFE [START_REF] Yin | A high-resolution 3D dynamic facial expression database[END_REF] and BP-4D-Spontaneous [START_REF] Zhang | BP4D-spontaneous: a high-resolution spontaneous 3D dynamic facial expression database[END_REF] datasets used a structured-light stereo sensor to capture textured 3D meshes of faces in controlled conditions and with various backgrounds. 2295 meshes were selected from these datasets and manually annotated with 66 landmarks and with self-occlusion information. Then, 16065 2D views were synthesized (seven views for each mesh) with yaw and pitch rotations ranging in the intervals [-45 • , +45 • ] and [-30 • , +30 • ], respectively. Additionally, there were 7,000 frames from the Multi-PIE (Gross et al, 1 Supplemental material for this paper can be found at https://team. inria.fr/perception/research/upa3dfa/ 2010) and 541 frames from the Time-Sliced datasets, respectively. Both these datasets contain RGB images gathered with multiple cameras from different viewpoints but with no 3D information. Therefore, a 3D face model is extracted for each image, using the model-based multi-view structure-from-motion technique of [START_REF] Jeni | Dense 3D face alignment from 2D video for real-time use[END_REF]. As above, each 3D face model was annotated with 66 landmarks and with self-occlusion information.

The following metrics were used in [START_REF] Jeni | The first 3D face alignment in the wild (3DFAW) challenge[END_REF]: the ground-truth error (GTE) and the cross-view ground-truth consistency error (CVGTCE), namely,

GTE(i) = 1 N N n=1 x n,i -xn,i /d i , (1)
CVGTCE(i) = 1 N N n=1 s i R i x n,i + t i -xn,i /d i , (2)
where x n,i denotes the n-th detected 3D landmark associated with test sample i, xn,i is the corresponding ground-truth 3D landmark, d i is the inter-ocular distance of the sample face i, N is the number of landmarks, and s i , R i , and t i are the scale factor, rotation matrix and translation vector associated with a rigid mapping that compensates the possible discrepancy between the set of detected landmarks and the set of ground-truth landmarks.

The Menpo challenge [START_REF] Deng | The Menpo benchmark for multi-pose 2D and 3D facial landmark localisation and tracking[END_REF] is based on a dataset of over 12000 face images. In order to obtain 2D and 3D ground-truth landmarks, an automatic annotation process is proposed, which fits a 3D face model to each 2D image. This fitting is carried out via non-linear minimization over the shape parameters (identity and expression), the rigid parameters (rotation and translation of the 3D model with respect to the camera), and the camera parameter (the scale of the weak-perspective model). The evaluation metric uses (1) with a different normalization factor, namely the size of the face bounding box.

The NoW benchmark is proposed in [START_REF] Sanyal | Learning to regress 3D face shape and expression from an image without 3D supervision[END_REF] for the task of 3D reconstruction from a single monocular image of a face. The associated dataset contains 2054 face images in frontal and profile views of 100 subjects and a 3D head scan for each subject. This dataset is similar in spirit to [START_REF] Bagdanov | The Florence 2D/3D hybrid face dataset[END_REF]. While the images contain four categories (neutral, expression, occlusion, and selfie) the 3D scans correspond to neutral faces. Therefore, the challenge for all categories is the reconstruction a neutral 3D face, which implies that non-neutral faces must undergo some form of disentanglement. Moreover, since the predicted 3D mesh and the ground-truth 3D scan lie in different coordinate systems, a rigid alignment must be performed prior to computing an evaluation metric. The authors provide an alignment procedure that minimizes a scan-to-mesh (or point-to-surface) distance over the alignment parameters (scale, rotation, and translation). This is a difficult alignment problem that necessitates to alternate between (i) selecting the closest points on the mesh and (ii) estimating the rigid parameters. Once an optimal alignment is found, the evaluation metric consists of the scan-to-mesh distance.

The evaluation metrics used in these benchmarks require ground-truth either of 3D landmarks or of 3D scans. Manual annotations of thousands of images cannot be achieved and automatic annotation must therefore be used. As outlined above, automatic annotation is based on complex non-linear minimization methods that are prone to errors and may not be reliable in the presence of profile views and of occlusions. Localization noise is inherent. Nevertheless, these evaluation metrics are limited in scope since they cannot distinguish between landmark localization noise (inlying data) and large localization errors (outlying data).

In contrast, the proposed methodology does not make use of ground-truth annotations. Robust rigid alignment (analyzed in detail below) is used to build frontal landmark models in a completely unsupervised way. A statistical characterization of each landmark is provided by measuring the discrepancy between the predicted landmark and the corresponding model landmark. Indeed, a confidence score is computed for each predicted landmark in order to assess its localization accuracy and to decide whether the landmark should be treated as an inlier (affected by detection noise) or as an outlier (the detection has failed). Our method may well be viewed as an analysis of performance of 3DFA algorithms, rather than a benchmark or a challenge. This is particularly useful whenever the output of 3DFA is used for facial expression recognition, for lip reading, for 3D pose estimation, etc. The proposed landmark analysis can also be applied to ground-truth landmarks in order to remove bad annotations, be them manual, semi-automatic or fully automatic.

A fundamental building block of the proposed methodology is the estimation of the rigid transformation (scale, rotation and translation) between two 3D point sets. We propose to perform this estimation in a robust way, where the term robustness refers to the capacity of a method to be unaffected by outliers. For that purpose, we cast the problem at hand in the framework of MLE, i.e. (7). In MLE, the choice of the likelihood function is crucial. We opt for two choices, namely a mixture model with two components, a Gaussian component and a uniform one (GUM) [START_REF] Banfield | Model-based gaussian and non-gaussian clustering[END_REF][START_REF] Zaharescu | Robust factorization methods using a Gaussian/uniform mixture model[END_REF][START_REF] Lathuilière | DeepGUM: Learning deep robust regression with a Gaussian-uniform mixture model[END_REF] and the generalized Student's t-distribution [START_REF] Mclachlan | Robust mixture modelling using the t distribution[END_REF][START_REF] Sun | Robust mixture clustering using Pearson type VII distribution[END_REF][START_REF] Forbes | A new family of multivariate heavy-tailed distributions with variable marginal amounts of tailweight: application to robust clustering[END_REF]. These two distributions behave quite differently. For each data point, GUM evaluates the posterior probability of being an inlier or an outlier, i.e. ( 16) and ( 18). Student evaluates a weight w associated with each data point. These weights play the role of precisions (the inverse of the variance): the higher the weights, the more reliable data. Each weight is treated as a random variable modeled with a gamma distribution, i.e. ( 24). Generalized Student belongs to the larger class of heavy-tailed distributions that are wellknown to be robust to outliers [START_REF] Mclachlan | Robust mixture modelling using the t distribution[END_REF]. To the best of our knowledge, this is the first time that a heavytailed distribution is used in conjunction with rigid alignment.

Both GUM and Student evaluate a posterior covariance matrix, i.e. ( 23) and ( 29), respectively. The evaluation of these covariance matrices is fundamental for taking into account heterogeneous landmark distributions, e.g. Fig. 1-right, and for as-sessing how much one can trust the results. This stays in strong contrast with the prevailing methods that are used in computer vision to estimate the rigid transformation between two point sets [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF][START_REF] Horn | Closed-form solution of absolute orientation using orthonormal matrices[END_REF][START_REF] Faugeras | The representation, recognition, and locating of 3-d objects[END_REF][START_REF] Arun | Least-squares fitting of two 3-D point sets[END_REF][START_REF] Umeyama | Least-squares estimation of transformation parameters between two point patterns[END_REF]. These methods assume an isotropic covariance matrix, i.e. spherical uncertainties, thus leading to closed-form solutions for estimating the rigid parameters. The fact that both GUM and Student estimate full covariance matrices, i.e. ellipsoidal uncertainties, introduces an additional complexity, namely a non-linear solver is required to estimate the rotation matrix (please consult Section 4.1). In [START_REF] Horaud | Rigid and articulated point registration with expectation conditional maximization[END_REF] this was addressed via convex relaxation. In this paper, we simplify the optimization problem by using quaternions and propose to use sequential quadratic programming [START_REF] Bonnans | Numerical optimization: theoretical and practical aspects[END_REF]) (please refer to Section 4.4 for an in-depth discussion).

The proposed confidence test detects anomalous outputs predicted by a 3DFA algorithm, be it based on a deep architecture or not. Alternatively, a number of methods were recently proposed to detect anomalous inputs, most notably outof-distribution (OOD) detection methods, e.g. [START_REF] Hendrycks | A baseline for detecting misclassified and outof-distribution examples in neural networks[END_REF][START_REF] Liang | Enhancing the reliability of out-of-distribution image detection in neural networks[END_REF][START_REF] Lee | Training confidence-calibrated classifiers for detecting out-of-distribution samples[END_REF][START_REF] Hendrycks | Deep anomaly detection with outlier exposure[END_REF]. These methods are concerned with classification networks that use a softmax layer to predict a probability distribution over the discrete set of possible labels (classes). Moreover these methods require annotated OOD examples. In contrast our proposed confidence test applies to regression networks and is unsupervised.

Problem Formulation

Let us consider the mapping between two sets of 3D facial landmarks, a predicted set, x 1:N = (x 1 . . . x N ) ∈ R 3×N , and a model set, y 1:N = (y 1 . . . y N ) ∈ R 3×N . The predicted set corresponds to a face with arbitrary and unknown pose, identity, expression and occlusion. Without loss of generality, the model set corresponds to a statistical frontal landmark model, namely each landmark n in this set is characterized by a mean vector y n and a covariance matrix C n ∈ R 3×3 , i.e. Section 6.2. In the general case, the unknown-to-frontal mapping is composed of a rigid transformation, i.e. head motion, and of a non-rigid facial deformation, up to an unknown error. Therefore, we introduce an additive residual r n ∈ R 3 , associated with the nth landmark, to account for the non-rigid component of the mapping and for various sources of errors. The mapping, from an unknown face pose to a frontal face pose can be modeled in the following way:

y n = sRx n + t + r n , ∀n ∈ {1, . . . , N}, (3)
where the rigid transformation is parameterized by a scale factor s ∈ R, a rotation matrix R ∈ S O(3) ⊂ R 3×3 and a translation vector t ∈ R 3 , while the non-rigid deformations and errors are represented by the residuals r 1:N . Let us assume that the optimal scale, rotation and translation are obtained via MLE, or equivalently, minimization of the negative log-likelihood pa-rameterized by θ:

{s , R , t , θ } = argmin s,R,t,θ N n=1 -log P(r n ; θ), (4)
followed by an estimation of the residuals:

r n = y n -s R x n -t , ∀n ∈ {1, . . . , N}. (5)
Note that the minimizer ( 4) must be immune to non-rigid deformations and to noise, which are jointly referred to as inliers, and at the same time it must be robust to large errors in landmark localization due to various perturbations, such as occlusions, motion blur, etc., which are referred to as outliers, i.e. Section 4. Therefore, we seek a robust rigid mapping technique that enables us to discriminate between inliers and outliers, namely:

       if r n C n ≤ 1 n = inlier otherwise n = outlier, ( 6
)
where C n = 9C n and • C is the Mahalanobis distance associated with covariance C. This guarantees with 99% confidence that if the nth landmark lies inside the ellipsoid defined by C n , then it is an inlier, and otherwise it is an outlier (it lies outside the ellipsoid), i.e. Section 6.3. The aptitude to discriminate between inlying and outlying landmarks is a core feature of the proposed unsupervised performance analysis.

Robust Rigid Mapping

We cast the problem of estimating the rigid mapping between two sets of landmarks into the framework of robust probability distribution functions. We assume that the residuals r 1:N are independent and identically distributed (i.i.d). Then, the problem of estimating the rigid transformation parameters can be solved via MLE or, equivalently, via negative log-likelihood minimization, namely min θ L(θ|(x, y) 1:N ), with:

L(θ|(x, y) 1:N ) = - N n=1 log P(r n ; θ), (7)
where P(r; θ) is the probability distribution function (pdf) of r parameterized by θ which is composed of s, R, t and of the pdf parameters that are specified below.

Gaussian Model

The simplest statistical model is to assume that the residuals follow a zero-centered Gaussian distribution with covariance matrix Σ ∈ R 3×3 , namely P(r; θ) = N(r; 0, Σ). By developing (7) and ignoring terms that do not depend on the model parameters, we obtain:

L(θ|(x, y) 1:N ) = 1 2 N n=1 y n -sRx n -t 2 Σ + log |Σ| , ( 8
)
where a 2 Σ = a Σ -1 a is the squared Mahalanobis norm of a ∈ R 3 . The minimization of (8) over t yields:

t = y -s R x, (9)
where the over-script indicates the optimal value of a parameter and with the notations:

x = N n=1 x n N , y = N n=1 y n N . ( 10
)
By substituting ( 9) into (8) and by using centered coordinates, i.e. x n = x n -x, y n = y n -y, we obtain:

L(θ|(x , y ) 1:N ) = 1 2 N n=1 y n -sRx n 2 Σ + log |Σ| . ( 11
)
Standard approaches to the minimization of ( 11) with respect to the rotation matrix assume an isotropic covariance, Σ = σI 3 . Indeed, the development of (11) yields

N n=1 sx n RΣ -1 R x n = sσ -1 N n=1
x n x n thus leading to closed-form solutions, e.g. [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF][START_REF] Horn | Closed-form solution of absolute orientation using orthonormal matrices[END_REF][START_REF] Faugeras | The representation, recognition, and locating of 3-d objects[END_REF][START_REF] Arun | Least-squares fitting of two 3-D point sets[END_REF][START_REF] Umeyama | Least-squares estimation of transformation parameters between two point patterns[END_REF] and Appendix Appendix A. Nevertheless, the isotropic-covariance assumption is barely valid in practice. In the case of a full covariance, the optimization becomes

R = argmin R 1 2 tr Σ -1 (s 2 RAR -2sRB) , (12)
where tr(•) is the trace operator and with the notations

A = N n=1 x n x n , B = N n=1
x n y n . A rotation matrix must satisfy RR = I 3 and |R| = +1. This yields a constrained non-linear optimization problem. An elegant formulation consists of parameterizing the rotation with a unit quaternion, thus reducing the number of parameters from nine to four, while the number of constraints is reduced from seven to one. Let R(q), where q is a unit quaternion (please consult Appendix Appendix A). Using this representation, the rotation is described by four parameters and the associated constrained optimization problem writes:

q = argmin q 1 2 tr Σ -1 (s 2 R(q)AR (q) -2sR(q)B) + λ q q -1 2 . ( 13
)
Similar to [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF] the optimal scale factor has a closedform expression:

s =       N n=1 y n Σ -1 y n N n=1 (R x n ) Σ -1 R x n       1/2 . ( 14
)
Finally, the optimal covariance is estimated with:

Σ = 1 N N n=1 (y n -s R x n )(y n -s R x n ) . ( 15
)
Once the rotation and scale are initialized using the method of [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF], alternating optimization can be used by iterating between ( 13), ( 14) and ( 15). We will refer to this method as generalized Horn.

Gaussian-uniform Mixture Model

Unfortunately, the above statistical model does not behave well in the presence of large residuals, or outliers. For the purpose of explicitly modeling inliers and outliers, a discrete random variable Z n is associated with each residual r n , and let z be a realization of Z. Now, r is drawn either from a zero-centered Gaussian distribution, as above, or from a multivariate uniform distribution:

P(r|Z = z) =        N(r; 0, Σ) if z = inlier U(r; 0, γ) if z = outlier, ( 16
)
where γ is the volume of the distribution. This yields a twocomponent mixture model, an inlier component with prior probability p, and an outlier component with prior probability 1p. This naturally leads to solving the problem via expectation-maximization (EM) which alternates betwen (i) evaluating the posterior probabilities of the residuals to be inliers or outliers, and (ii) minimizing the expected complete-data negative log-likelihood, E Z [-log P(r 1:N , Z 1:N |r 1:N ; θ)], where the expectation is taken over the realizations of Z, and where the parameter vector is θ = {s, R, p, Σ}. 2 This yields the minimization of:

E(θ|(x , y ) 1:N ) = 1 2 N n=1 α n y n -sRx n 2 Σ + log |Σ| (17)
where the posterior probability to be an inlier α n = P(Z = inlier|r n ), is :

α n = p N(r n ; 0, Σ) p N(r n ; 0, Σ) + (1 -p) γ -1 , (18)
and the posterior probability to be an outlier is P(Z = outlier|r n ) = 1 -α n . The presence of α 1:N in ( 17) replaces (10) with:

x = N n=1 α n x n N n=1 α n , y = N n=1 α n y n N n=1 α n , (19)
as well as A and B from (13) with

A = N n=1 α n x n x n , B = N n=1 α n x n y n . (20)
Hence, (13) can be used to estimate the optimal rotation. Moreover, ( 14) is replaced with :

s =       N n=1 α n y n Σ -1 y n N n=1 α n (R x n ) Σ -1 R x n       1/2 . ( 21
)
The prior probability p and covariance matrix Σ are estimated with:

p = 1 N N n=1 α n , (22)
Σ = N n=1 α n (y n -s R x n )(y n -s R x n ) N n=1 α n . ( 23
)
We refer to this model as the Gaussian-uniform mixture (GUM) and the associated EM is summarized in Algorithm 1.

Data: Centered point coordinates, i.e. (10). Normalization parameter γ ; Initialization of θ old : Use the closed-form solution [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF] to evaluate s old and R old and then use these parameter values to evaluate Σ old using (15) and set p old = 0.8; while θ newθ old > do E-step: Evaluate the posteriors α 1:N using (18) with θ old ; Update the centered coordinates using (19) ; M-scale-step: Evaluate s new using (21); M-rotation-step: Estimate R new via constrained non-linear optimization of ( 13) using ( 20) ; M-covariance-step: Evaluate Σ new using (23); M-prior-step: Evaluate p new using ( 22); θ old ← θ new ; end Optimal translation: Evaluate the translation vector using (9); Result: Estimated scale s , rotation R , translation t , prior p , covariance Σ , and posterior probabilities of landmarks α 1:N . Algorithm 1: GUM-EM for robust estimation of the rigid transformation between two point sets.

Generalized Student Model

Another way to enforce robustness is to use the generalized Student's t-distribution, also known as the Pearson type VII distribution [START_REF] Sun | Robust mixture clustering using Pearson type VII distribution[END_REF]:

P(r; Σ,µ, ν) = ∞ 0 N(r; 0, w -1 Σ)G(w, µ, ν)dw = Γ(µ + 3 2 ) |Σ| 1 2 Γ(µ)(2πν) 3 2        1 + r 2 M 2ν        -(µ+ 3 2 ) , ( 24
)
where µ and ν are the parameters of the prior gamma distribution of w and Γ(•) is the gamma function. The distribution (24) differs from the standard Student's t-distribution in that the weight variable W, or the precision, is drawn from a gamma distribution with parameters µ and ν, instead of ν/2 and ν/2. Notice that in (24) ν and Σ appear only through their product, which means that an additional constraint is required to make the parameterization unique. One possibility is to constrain the determinant of the covariance, e.g. |Σ| = 1, which is equivalent to have an unconstrained Σ with ν = 1. Unconstrained parameters are easier to deal with in inference algorithms. Therefore, we will rather assume without loss of generality that ν = 1.

Notice that the posterior distribution of w is also a gamma distribution, namely the posterior gamma distribution:

P(w n |r n ; Σ, µ, ν) = N(r n ; 0, w -1 n Σ)G(w n , µ, ν) = G(w n ; a, b n ), (25)
with parameters:

a = µ + 3 2 , b n = 1 + 1 2 r n 2 Σ . (26)
Data: Centered point coordinates, i.e. (10). ; Initialization of θ old : Use the closed-form solution [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF] to evaluate s old and R old ; evaluate Σ old using (15). Provide µ old ; while θ newθ old > do E-step: evaluate a new and b new 1:N using (26) with θ old , then evaluate w new 1:N using (27) ; Update the centered coordinates using ( 19), where α 1:N are replaced with w 1:N ; M-scale-step: Evaluate s new using (21); M-rotation-step: Estimate R new with ( 13), (20) ; M-covariance-step: Evaluate Σ new using (29) ; M-mu-step: Evaluate µ new using (30) ; θ old ← θ new ; end Optimal translation: Evaluate the translation vector using (9); Result: Estimated scale s , rotation R , translation t , covariance Σ , and landmark weights w 1:N . Algorithm 2: The GStudent-EM for robust estimation of the rigid transformation between two point sets.

The posterior mean of the weight variable is:

w n = E[w n |r n ] = a b n . (27)
As with the Gaussian-uniform model, we need to minimize the expected complete-data negative log-likelihood, E W [-log P(r 1:N , W 1:N |r 1:N ; θ)] and in this case the parameter vector is θ = {s, R, Σ, µ} since we set ν = 1. This yields the minimization of:

Q(θ|(x , y ) 1:N ) = 1 2 N n=1 w n y n -sRx n 2 Σ + log |Σ| , (28)
thus replacing α 1:N with w 1:N in ( 19) and ( 20) to estimate the optimal rotation (13) and scale (21). The covariance matrix is estimated with:

Σ = 1 N N n=1 w n (y n -sRx n )(y n -sRx n ) (29)
The parameter µ is updated by solving the following equation, where Ψ(•) is the digamma function:

µ = Ψ -1        Ψ(a) - 1 n N n=1 log b n        . ( 30
)
We refer to this model as the generalized Student (GStudent) and the associated EM algorithm is summarized in Algorithm 2.

Algorithm Implementation and Analysis

Algorithm 1 and Algorithm 2 are expectation maximization (EM) procedures and it is well known that they have good convergence properties. One should notice that all the computations inside these algorithms are in closed-form, with the notable exception of the estimation of the rotation matrix. The latter is parameterized with a unit quaternion which may be estimated via nonlinear constrained optimization. The unitquaternion parameterization of rotations, i.e. Appendix Appendix A, has several advantages: (i) the number of parameters to be estimated is reduced from nine to four, (ii) the number of nonlinear constraints is reduced from seven constraints (six quadratic constraints, i.e. R R = I, and one quartic constraint, i.e. |R| = 1) to one quadratic constraint (q q = 1), (iii) the initialization is performed with the closed-form solution of [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF] that uses a unit quaternion as well.

In practice, the constrained nonlinear optimization problem ( 13) is solved using the sequential quadratic programming method [START_REF] Bonnans | Numerical optimization: theoretical and practical aspects[END_REF], more precisely a sequential least squares programming (SLSQP) solver3 is used in combination with a root-finding software package [START_REF] Kraft | A software package for sequential quadratic programming[END_REF]. The SLSQP minimizer found at the previous EM iteration is used as an initial estimate at the current EM iteration. The closed-form method of [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF] (please consult Appendix Appendix A) is used to initialize the unit-quaternion parameters at the start of the EM algorithm.

Other closed-form methods commonly used in computer vision, e.g. [START_REF] Horn | Closed-form solution of absolute orientation using orthonormal matrices[END_REF][START_REF] Arun | Least-squares fitting of two 3-D point sets[END_REF], perform singular value decomposition to extract an orthogonal matrix from the measurement matrix, but without the guarantee that the estimated matrix is a rotation, i.e. its determinant must be equal to +1. Appendix Appendix A summarizes the unit-quaternion closed-form method, which is based on estimating the smallest eigenvalue and eigenvector pair of a 4×4 semi-definite positive symmetric matrix -a well known mathematical problem yielding a straightforward numerical solver.

Analyzing the Robustness of Rigid Mapping

In order to quantify the performance of the proposed robust rigid-mapping algorithms, we devised an experimental protocol on the following grounds. Let x 1:N be a set of landmarks associated with the frontal view of a face. The set y m 1:N is generated with:

y m n (b) = s m R m x n + t m + r m n (b), ∀n ∈ {1, . . . , N}, (31)
where b > 0 is a scalar that controls the level of noise and m is the trial index. As described in detail below, the noise level, b can be the variance of Gaussian isotropic noise, the total variance of Gaussian anisotropic noise, or the volume of uniformly distributed noise. The landmark coordinates are normalized such that ∀n, x n ∈ [0, 1] 3 . For each noise level, we randomly generate M trials, namely M rigid mappings and M sets of N residuals r 1:N = {r n } n=N n=1 . For each trial m we estimate the rigid mapping parameters, s m , R m , t m , and we measure the root mean square error (RMSE) between these estimated parameters and the ground-truth parameters, sm , Rm , tm , namely:

E s = 1/M M m=1 |s m -sm | 2 1/2 , ( 32
) E t = 1/M M m=1 t m - tm 2 1/2 , ( 33
) E R = 1/M M m=1 R m - Rm 2 1/2 , (34)
The ground-truth rigid-mapping parameters are generated in the following way. For each trial m, the scale and the translation vector are generated from uniform distributions, namely s m ∼ U(0.5, 2) and t m ∼ U(0.5, 5) 3 . The rotation matrix is parameterized by the pan, tilt and yaw angles, namely:

R = R γ R φ R ψ =           cos γ -sin γ 0 sin γ cos γ 0 0 0 1                     cos φ 0 sin φ 0 1 0 -sin φ 0 cos φ                     1 0 0 0 cos ψ -sin ψ 0 sin ψ cos ψ          
A rotation matrix is obtained by randomly generating the pan, tilt and yaw angles, γ m , φ m , ψ m , from a uniform distribution, U(-90 • , +90 • ).

In order to generate residuals, r 1:N , we simulate three types of noise:

• Isotropic Gaussian noise: r ∼ N(0, σI);

• Anisotropic Gaussian noise: r ∼ N(0, Σ), and

• Uniform noise: r ∼ U(-a/2, a/2) 3 .
In the case of anisotropic noise, a covariance matrix must be randomly generated for each trial. This is done in the following way. Let Σ = QΛQ , with Q ∈ O(3) (an orthogonal matrix) and with Λ = Diag (λ 1 , λ 2 , λ 3 ), where the eigenvalues correspond to the variances along the eigenvectors -the directions of maximum variance. Let λ = λ 1 + λ 2 + λ 3 denote the total variance. A sample covariance matrix Σ is simulated by randomly generating an orthogonal matrix Q and by randomly generating the three eigenvalues from a uniform distribution, U(0, 1).

We tested the following rigid mapping models and associated algorithms:

• Horn: Gaussian distribution with isotropic covariance, [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF] and Appendix Appendix A;

• Gen-Horn: Gaussian distribution with anisotropic covariance, Section 4.1;

• GUM-EM: Gaussian-uniform mixture distribution, Algorithm 1, and

• GStudent-EM: Generalized Student's t-distribution, Algo- rithm 2.
The experiments were conducted in the following way. For each noise level, we simulated M = 500 trials for which we computed the RMSEs, namely eqs. ( 32), (33), and (34). For each trial m we split the landmarks into an inlier set and an outlier set and the N = 68 landmarks are randomly assigned to one of these sets. The first experiment determines the percentage of outliers that can be handled by the robust algorithms, Figure 2. For this purpose, the percentage of outliers is increased from 10% to 60%. The inlier noise is drawn from an anisotropic Gaussian distribution with a total variance λ = 0.0025. The outlier noise is drawn from a uniform distribution with amplitude a = 1.5 (remember that the landmark coordinates are normalized to lie in the interval [0, 1]). The cuves plotted in Figure 2 show that the RMSE associated with non robust methods, i.e. Horn and Gen-Horn increase monotonically. On the contrary, the robust algorithms, GUM-EM and GStudent-EM, have a radically different behavior. After a short increase, the RMSE remains constant, and then it increases again.

In the other experiments, the number of inliers was set to be equal to the number of outliers and we experimented with the three noise types already mentioned. Figure 3 shows the RMSEs when inlier noise is drawn from an isotropic Gaussian distribution with σ = 0.0025, while outlier noise is drawn from a uniform distribution whose volume is increased from a = 0.2 to a = 1.0. Similarly, Figure 5 shows the RMSEs for the case when inlier noise is drawn from an anisotropic Gaussian distribution with total variance λ = 0.0025, while outlier noise is drawn from a uniform distribution whose volume is increased from a = 0.2 to a = 1.0. Finally, Figure 5 shows the RM-SEs when inlier noise is drawn from an anisotropic Gaussian distribution with total variance λ = 0.0025, while outlier noise is drawn from an anisotropic Gaussian distribution with total variance varying from λ = 0.2 to λ = 1.0.

These experiments clearly show that the two classes of methods (non-robust and robust) behave differently. The performance of non-robust rigid mapping decreases monotonically in the presence of outliers with increasing levels of noise. The robust methods can deal with up to 50% of outliers affected by a substantial noise level (1.5

times the size of the image).

There is no evidence that the Gen-Horn algorithm performs better than the standard Horn algorithm. Nevertheless, Gen-Horn provides interesting information about the 3D structure of the estimated anisotropic covariance. The GUM-EM algorithm performs slightly better than the GStudent-EM algorithm, in particular in the presence of outliers drawn form a uniform distribution.

Measuring the Performance of 3D Face Alignment

In this section we describe an unsupervised methodology for quantitatively assessing the performance of 3DFA algorithms. The idea of the proposed benchmarking is to apply 3DFA to a dataset of face images in order to extract 3D landmarks, to robustly estimate the rigid transformation that maps these facial landmarks into a 3D landmark model, and to analyze the discrepancy between the extracted 3D landmarks and the model. Based on a confidence score, it is then possible to decide whether a landmark is correctly localized or not. This allows to assess the overall performance of a 3DFA algorithm as well as its behavior with respect to various perturbations, such as occlusions or motion blur.

Neutral Frontal Landmark Model

We start by computing a neutral frontal landmark model y 1:N in the following way. For this purpose, we use a dataset D 1 of K images of neutral faces (frontal viewing, no expression and no interfering object causing occlusion) and we extract N landmarks from each one of these K faces, {y 1:N,k } k=K k=1 . Then we use the landmark coordinates to compute the directions of maximum variance (or the principal components) of each face. By aligning these directions over the dataset, we compute a mean for each landmark, namely

y n = 1/K K k=1 y n,k , ∀n, 1 ≤ n ≤ N.
(35)

Statistical Frontal Landmark Model

We now explain how a statistical frontal landmark model is built, namely { p 1:N , C 1:N }, where p 1:N is the set of N means and C 1:N is the set of N covariance matrices associated with the statistical frontal landmark model. For this purpose, we use another dataset D 2 that contains L images of faces with the following characteristics: arbitrary poses, arbitrary expressions, both speaking and silent faces, but with no external sources of perturbation such as the presence of interfering object that may cause occlusions. We extract 3D landmarks from these L images using a 3DFA algorithm, namely {x 1:N,l } l=L l=1 , and we use either GUM-EM (Algorithm 1) or GStudent-EM (Algorithm 2) to robustly estimate the rigid transformations between each landmark-set l and the the neutral frontal landmark-set x 1:N,l and y 1:N . Based on this, we obtain L rigid-mapping parameters (one for each l): L scale factors, L rotations and L translations:

{s Alg l , R Alg l , t Alg l } L l=1
, where the over-script Alg denotes a robust algorithm, namely either GUM-EM or GStudent-EM. We remind that both algorithms provide a figure of merit for each landmark: posterior probabilities {α n,l } n=N n=1 in the case of GUM-EM, i.e. (18), and precision means {w n,l } n=N n=1 in the case of i.e. (27). Applying one of these robust rigid-alignment methods provides frontal landmarks, { x Alg 1:N,l } l=L l=1 , namely:

x

Alg n,l = s Alg l R Alg l x n,l + t Alg l . (36)
There are two different expressions for the posterior means and posterior covariances for GUM-EM and for GStudent-EM, respectively:

p GUM n = L l=1 α n,l x GUM n,l L l=1 α n,l , (37)
C GUM n = L l=1 α n,l ( x GUM n,l -p GUM n )( x GUM n,l -p GUM n ) L l=1 α n,l , (38)
and

p GSt n = L l=1 w n,l x GSt n,l L l=1 w n,l , (39)
C GSt n = L l=1 w n,l ( x GSt n,l -p GSt n )( x GSt n,l -p GSt n ) L l=1 w n,l . (40)
Notice that ( 39) and ( 40) compute a mean and a covariance for landmark n over the entire dataset. Hence, and unlike in (29), the covariance should be normalized with the sum of the weights.

Unsupervised Confidence Test

We now develop an unsupervised (statistical) confidence test for assessing whether the accuracy of a landmark, i.e. its 3D coordinates, is within (inlier) or outside (outlier) an expected range [START_REF] Savage | The Foundations of Statistics[END_REF]. Let us drop the algorithm over-script and let C n = Q n Λ n Q n be the eigen factorization of C n , where Q n is an orthonormal matrix and Λ n is a diagonal matrix containing the eigenvalues. We can now project each landmark (n, l) onto the space spanned by the three eigenvectors of this matrix:

z n,l = Q n ( x n,l -p n ) (41)
Landmark (n, l) is an inlier with 99% confidence if z n,l lies inside the ellipsoid whose half-axes are three times the standard deviations, or 3

λ 1 n , 3 λ 2 n , 3 λ 3 n , where {λ 1 n , λ 2 n , λ 3 n } are the eigenvalues of C n , or z n,l Λ -1 n z n,l ≤ 1 ( 42
)
where Λ n = 9Λ n . Combining ( 41) and (42), yields ( x n,l -

p n ) Q n Λ -1 n Q n ( x n,l -p n ) ≤ 1. With the notation C n = Q n Λ n Q n . (43)
The 99% confidence test writes:

       if x n,m -p n C n ≤ 1 (n, m) = inlier otherwise (n, m) = outlier (44)
Based on this confidence test, we can now build a confidencetest accuracy (the higher the better) associated with a sample face m, namely:

u(m) = 1 N N n=1 I( x n,l -p n C n ≤ 1), (45)
where I(•) denotes the indicator function. Notice that (45) corresponds to the percentage of inliers, i.e. landmarks that, once scaled, rotated and translated, lie inside the confidence volume. Therefore, (45) can be used to assess whether the pose has been correctly estimated, namely u ≤ 50%, or not (please consult Section 5). For a test dataset D 3 composed of M samples, one can then compute the mean confidence test accuracy (the higher the better):

U = 1 M M m=1 u(m) (46)

Correlation with Supervised Metrics

In general, datasets of faces come with their ground-truth annotations, and we denote with x1:N,1:M the set of ground-truth landmarks associated with the dataset D 2 . We modify ( 1) to be able to build a metric that counts the proportion of inliers, namely the ground-truth accuracy (the higher the better):

s(m) = 1 N N n=1 I( x n,m -xn,m /d m ≤ ε), (47)
where ε is a user-defined threshold that corresponds to the quality of the ground-truth landmarks. Based on this we can compute the mean ground-truth accuracy (S):

S = 1 M M m=1 s(m) (48)
Finally, another interesting metric is the correlation coefficient between the above unsupervised and supervised metrics:

Cor = M m=1 (u(m) -U)(s(m) -S) M m=1 (u(m) -U) 2 M m=1 (s(m) -S) 2 1/2
(49)

Experimental Results

Neutral Frontal Landmark Model

The neutral frontal landmark model was trained in-the-wild by harvesting web images and using a face detector and a headpose estimator in order to select frontal faces. These images were visually inspected to guarantee shape and aspect variabilities as well as neutral facial expressions. This process yields a dataset D 1 composed of 1, 000 images. We used the 3DFA method of [START_REF] Feng | Joint 3D face reconstruction and dense alignment with position map regression network[END_REF] to extract landmarks from each face in the dataset. Next, we aligned them (please consult Section 6.1) and computed the landmark means using (35). Figure 6 shows a few examples of images from this dataset as well as the detected landmarks. Figure 7 show the neutral frontal landmark model thus obtained.

Statistical Frontal Landmark Model

The statistical frontal landmark model was trained from the YawDD dataset [START_REF] Abtahi | Yawdd: A yawning detection dataset[END_REF]. This dataset contains 322 videos which is equivalent to approximatively 300, 000 images. The face images in this dataset have large variabilities in terms of face shapes, face aspects, head poses and facial expressions. All the images were processed with no human intervention, namely: face detection, 3D face alignment, and robust rigid alignment with the neutral face landmarks just described. This yields the statistical face landmark model described in Section 6.2. For that purpose we used two 3DFA methods and the two robust alignment algorithms described in this paper. Hence, there are four possible 3DFA and robust alignment combinations that we used to train four different models:

• 3DFA1/GUM-EM: [START_REF] Bulat | Two-stage convolutional part heatmap regression for the 1st 3D face alignment in the wild (3DFAW) challenge[END_REF] and GUM-EM (Algorithm 1),

• 3DFA2/GUM-EM: [START_REF] Feng | Joint 3D face reconstruction and dense alignment with position map regression network[END_REF] and GUM-EM (Algorithm 1),

• 3DFA1/GStudent-EM: [START_REF] Bulat | Two-stage convolutional part heatmap regression for the 1st 3D face alignment in the wild (3DFAW) challenge[END_REF] and GStudent-EM (Algorithm 2), and

• 3DFA2/GStudent-EM: [START_REF] Feng | Joint 3D face reconstruction and dense alignment with position map regression network[END_REF] and GStudent-EM (Algorithm 2).

Figure 8 shows the statistical frontal landmark models obtained with these four combination. In this figure, the dots correspond to the posterior means, i.e. ( 37) and ( 39), while the ellipses correspond to image projections of the ellipsoids defined by (43).

Performance Evaluation of 3D Face Alignment

Once the neutral-frontal and statistical-frontal models are computed using datasets D 1 and D 2 , respectively, we use a third dataset, D 3 , to empirically assess the performance of four 3DFA algorithms using the unsupervised confidence test introduced in Section 6.3:

• 3DFA1 (Bulat and Tzimiropoulos, 2016) (ECCVW'16) and [START_REF] Bulat | How far are we from solving the 2D & 3D face alignment problem? (and a dataset of 230,000 3D facial landmarks)[END_REF]) (ICCV'17);

• 3DFA2 [START_REF] Feng | Joint 3D face reconstruction and dense alignment with position map regression network[END_REF]) (ECCV'18);

• 3DFA3 [START_REF] Zhu | Face alignment across large poses: a 3D solution[END_REF]) (CVPR'18) and [START_REF] Zhu | Face alignment in full pose range: A 3d total solution[END_REF] (PAMI'19), and

• 3DFA4 (Tu et al, 2020) (IEEE TMM'20).

For this purpose we used AFLW2000-3D [START_REF] Zhu | Face alignment across large poses: a 3D solution[END_REF] Each image in the AFLW2000-3D dataset is annotated with 68 3D landmarks. This semi-automatic annotation is performed by fitting a 3D deformable model to a dataset of 2D face images, e.g. [START_REF] Ghiasi | Occlusion Coherence: Localizing Occluded Faces with a Hierarchical Deformable Part Model[END_REF]. Nevertheless, as noted in [START_REF] Bulat | How far are we from solving the 2D & 3D face alignment problem? (and a dataset of 230,000 3D facial landmarks)[END_REF], many of the annotated landmarks in this dataset have large localization errors, especially in the case of profile views. Hence, performance evaluation based on supervised metrics are prone to errors. In [START_REF] Bulat | How far are we from solving the 2D & 3D face alignment problem? (and a dataset of 230,000 3D facial landmarks)[END_REF] it is visually shown that in these extreme poses their 3DFA method yields more precise landmark localization than the automatically annotated ones. Based on these observations, we applied our unsupervised performance analysis to the annotated landmarks as well, yielding the following combinations: • GT/GUM-EM: Ground-truth landmarks provided by [START_REF] Zhu | Face alignment across large poses: a 3D solution[END_REF] and GUM-EM, and

• GT/GStudent-EM: Ground-truth landmarks provided by [START_REF] Zhu | Face alignment across large poses: a 3D solution[END_REF] and GStudent-EM.

The results based on computing the mean confidence-test accuracies, i.e. ( 46) are summarized in Table 1. We remind that we used different datasets for training the neutral and statistical face landmark models, i.e. D 1 and D 2 , and for assessing the performance of the various combinations of 3DFA methods and robust-rigid mappings, i.e. D 3 . The means, evaluated over the confidence-test scores obtained with the annotated (groundtruth) landmarks (the last two rows of Table 1), are equal to 0.70 and to 0.65, respectively, which seems to confirm that the ground-truth landmark locations in the AFLW2000-3D dataset contain a substantial amount of errors and that, overall, both 3DFA methods that we analyzed, [START_REF] Bulat | Two-stage convolutional part heatmap regression for the 1st 3D face alignment in the wild (3DFAW) challenge[END_REF] and [START_REF] Feng | Joint 3D face reconstruction and dense alignment with position map regression network[END_REF], predict landmark locations that are more accurate than the ground-truth locations themselves.

We now compute correlation coefficients, i.e. (49), between the unsupervised and supervised metrics, i.e. (45), and (47). The results are reported in Table 2. Notice however that the supervised scores depend on the choice of the parameter ε. As done in [START_REF] Bulat | How far are we from solving the 2D & 3D face alignment problem? (and a dataset of 230,000 3D facial landmarks)[END_REF], this parameter was adjusted to eliminate samples yielding a very low score. With ε = 0.1 in (48) we obtained the following scores:

• 3DFA1: S = 0.57,

• 3DFA2: S = 0.60,

• 3DFA3: S = 0.49,

• 3DFA4: S = 0.09.

These scores are comparable with the scores reported in [START_REF] Bulat | How far are we from solving the 2D & 3D face alignment problem? (and a dataset of 230,000 3D facial landmarks)[END_REF] which uses a different normalization parameter. Notice that the scores obtained with the proposed confidence test, i.e. Table 1, are higher than these scores.

In the light of these results, we attempted to analyse the effect of eliminating inaccurate ground-truth landmark annotations from the benchmark just described. Let us define

M τ = {m | û(m) ≥ τ}, M τ ⊂ D 3 , (50)
where û(m) denotes the value of the unsupervised score ( 45) associated with the ground-truth landmarks of face sample m.

We see that when τ increases, the accuracy of the ground-truth landmark annotations contained in the subset M τ increases as well, at the price of drastically decreasing the number of inlying samples, which in turn lowers down the statistical significance of the resulting scores. The correlation coefficient is then computed with the following formula:

Cor(τ) = m∈M τ (u(m) -U)(s(m) -S) m∈M τ (u(m) -U) 2 m∈M τ (s(m) -S) 2 1/2
(51)

Figure 9-left shows the correlation as a function of τ where the 3DFA1/GStudent-EM method was used to build the confidence test, i.e. Figure 8-c, while Figure 9-right shows the corresponding p-value (the smaller the better) with a significance level of 0.005: a very small p-value is an indicator of the statistical significance of the correlation measure. The red dots in these plots correspond to a p-value not satisfying the significance level just mentioned.

In the light of these experiments, we conclude that the proposed methodology for assessing the performance of 3DFA methods is not biased by the quality of landmark annotation, whether the latter is automatic or human-assisted. The experiments suggest that the proposed unsupervised methodology could be used (i) to assess the quality of landmark annotation itself and (ii) to remove badly annotated landmarks.

We now illustrate the proposed performance analysis method with a few examples from the AFLW2000-3D. In all these examples, the statistical face model is trained with 3DFA1/GStudent-EM and GStudent-EM is used for testing. Figure 10 shows results obtained with 3DFA1 [START_REF] Bulat | Two-stage convolutional part heatmap regression for the 1st 3D face alignment in the wild (3DFAW) challenge[END_REF], Figure 11 shows results obtained with 3DFA2, Figure 12 shows results obtained with 3DFA3, and Figure 13 shows results obtained with 3DFA4.

We now show results obtained by applying GStudent-EM to the ground-truth landmarks associated with the AFLW2000-3D dataset [START_REF] Zhu | Face alignment across large poses: a 3D solution[END_REF], or GT/GStudent-EM. Some best-score examples are shown in Figure 14 and some worse-score results are shown in Figure 15. Both the results reported in Table 1 and these examples suggest that the ground-truth annotations should be handled with caution.

We also gathered a dataset of 30 animal faces in order to analyse the performance of 3DFA with non-human faces. We processed these 30 images with the four 3DFA methods, as obove, but only [START_REF] Bulat | Two-stage convolutional part heatmap regression for the 1st 3D face alignment in the wild (3DFAW) challenge[END_REF] yielded exploitable results. As with human faces, for each one of these animal images we computed the percentage of inliers, i.e. (45). Figure 16 shows the six best and the six worst scores. The high scores observed on some of these animal faces correspond to failures of the proposed confidence test. In fact, the 3DFA method itself failed in these cases because it predicted a landmark pattern that corresponds to a human face. These results should, however, be interpreted with caution, because neither the 3DFA method of [START_REF] Bulat | Two-stage convolutional part heatmap regression for the 1st 3D face alignment in the wild (3DFAW) challenge[END_REF] nor the proposed methodology were trained with animal faces.

Discussion and Conclusions

Landmarks predicted by face alignment methods, whether 2D or 3D, are inherently affected by noise and they contain outliers. We propose an unsupervised methodology to characterize the performance of 3DFA. First, we learn a statistical frontal landmark model, namely posterior mean and covariance for each 3D landmark. This pose-invariant model is materialized by an ellipsoidal volume of confidence that encapsulates variabilities due to face identity, face expression, occlusion, and detection noise. Second, the landmarks predicted by a 3DFA method are rigidly and robustly mapped onto this frontal model in order to be able to compare their locations with the model locations. A landmark that falls inside its corresponding ellipsoidal volume is labeled as inlier with 99% confidence, while a landmark that lies outside this volume is labeled as outlier. The ability (i) to separate pose from expression and from identity in a robust manner and (ii) to discriminate between inlying and outlying landmarks stands in contrast with existing evaluation metrics that compute the mean distance between ground-truth landmark locations and predicted locations.

We note that none of the 3DFA evaluation metrics proposed so far exploit the concept of bringing all the facial landmarks into a canonical frontal pose. This is performed using a rigid alignment that is embedded into a maximum-likelihood estimator that uses a robust probability distribution function. This may well be viewed either as a robust pose estimator or as a mechanism that yields an expression-and identity-preserving frontal landmark representations. In turn, this enables temporal analysis of facial expressions and of lip movements.

When applied to the AFLW2000-3D dataset, the proposed analysis reveals that ground-truth landmark locations, provided by this annotated set of faces, contain 0.67 inliers, on an average (the last two rows of Table 1), which is less than the percentage of inliers associated with the three best-performing 3DFA methods that we analyzed. This result confirms the conclusions of [START_REF] Bulat | How far are we from solving the 2D & 3D face alignment problem? (and a dataset of 230,000 3D facial landmarks)[END_REF]) that these annotations contain many large errors. To better understand these results, we also computed the correlation between the proposed unsupervised metric and the supervised metric. The correlation coefficients are in between 0.12 and 0.33 for the best-performing methods, as reported in Table 2. Interestingly, these correlation coefficients increase monotonically as the bad annotations are eliminated, which is illustrated in Figure 9. We conjecture that the proposed methodology can be used to eliminate annotation errors from a dataset of faces. Alternatively, it is also possible to automatically annotate a dataset using existing 3DFA methods and selecting the best localization for each landmark.

Because the proposed methodological pipeline makes use of 3DFA and of robust rigid mapping both for training and for testing, one may argue that the associated performance metric is biased. We empirically showed that the analysis is agnostic to various combinations of methods used for training the model and for the tests themselves. Moreover, we showed that the method could also be used to assess the quality of facial landmark annotations, in particular those annotations that are obtained automatically based on 3D deformable-model fitting. We therefore conclude that the interest of the proposed methodology is twofold, namely (i) to assess the accuracy, the repeatability and the reliability of the predictions obtained with 3DFA algorithms, and (ii) to evaluate the quality of the landmarks predicted from a test face.

Table 1: Performance analysis based on the proposed unsupervised metrics. The numbers correspond to the proportion of inliers (the higher the better) computed using ( 45) and ( 46) over a dataset that contains 2,000 face images and 68 landmarks per face.

Alignment method

Statistical face models trained with D 1 and D 2 datasets: using dataset D

  (a) RMSE error in rotation. (b) RMSE error in scale. (c) RMSE error in translation.

Fig. 2 :

 2 Fig.2: RMSE error as a function of the percentage of outliers: inliers are affected by anisotropic Gaussian noise with total variance λ = 0.0025, while the percentage of outliers, affected by uniform noise with amplitude a = 1.5, increases from 0% to 60%.

Fig. 3 :

 3 Fig. 3: RMSE error as a function of uniform noise affecting a fixed number of outliers: inliers (50%) are affected by isotropic Gaussian noise with variance σ = 0.0025, while outliers (50%) are affected by uniform noise of increasing amplitude a ∈ [0.2, 1.0].

Fig. 4 :

 4 Fig. 4: RMSE error as a function of uniform noise affecting a fixed number of outliers: inliers (50%) are affected by anisotropic Gaussian noise with total variance λ = 0.0025, while outliers (50%) are affected by uniform noise of increasing amplitude a ∈ [0.2, 1.0].

  (a) RMSE error in rotation. (b) RMSE error in scale. (c) RMSE error in translation.

Fig. 5 :

 5 Fig. 5: RMSE error as a function of anisotropic Gaussian noise affecting a fixed number of outliers: inliers (50%) are affected by anisotropic Gaussian noise with total variance λ = 0.0025, while outliers (50%) are affected by anisotropic Gaussian noise with total variance λ ∈ [0.2, 1.0].

  as a test dataset, consisting of 2000 images with large-pose variations. More precisely, the yaw angles (vertical axis of rotation) in the following intervals [0 • , ±30 • ] for 1306 faces, in the interval [±30 • , ±60 • ] for 462 faces and in the interval [±60 • , ±90 • ] for 232 faces. The dataset contains a large variety of facial shapes, facial expressions and illuminations conditions. Moreover, there are many faces with partial occlusions caused by the presence of hair, hands, handheld objects, glasses, etc. Notice that large poses induce partial occlusions as well.

Fig. 6 :

 6 Fig. 6: Examples of faces and corresponding landmarks used to compute a neutral frontal landmark model.

Fig. 7 :

 7 Fig. 7: A 3D view of the neutral frontal landmark model.

Fig. 8 :

 8 Fig. 8: Statistical frontal landmark models obtained with two 3DFA methods and with the proposed robust rigid-mapping algorithms.

Fig. 9 :

 9 Fig.9: Correlation between the unsupervised and supervised metrics (left) and corresponding p-value (right) as a function of the accuracy of the ground-truth landmark annotations. The red dots correspond to a p-value not satisfying a significance level, which is set to 0.005 in these experiments.

Fig. 10 :

 10 Fig. 10: A few examples obtained with 3DFA1 (Bulat and Tzimiropoulos, 2016) and with GStudent-EM.

Fig. 11 :

 11 Fig.11:A few examples obtained with 3DFA2[START_REF] Feng | Joint 3D face reconstruction and dense alignment with position map regression network[END_REF] and with GStudent-EM.

Fig. 12 :

 12 Fig. 12: A few examples obtained with 3DFA3 (Zhu et al, 2016) and with GStudent-EM.

Fig. 14 :

 14 Fig. 14: Some examples of the best scores obtained with GStudent-EM and the ground-truth landmarks.

Fig. 15 :

 15 Fig. 15: Some examples of the worse scores obtained with GStudent-EM and the ground-truth landmarks.

Fig. 16 :

 16 Fig. 16: Percentage of inliers obtained with 12 animal faces using 3DFA1 (Bulat and Tzimiropoulos, 2016) (from top to bottom and left to right: 0.91, 0.88, 0.88. 0.87, 0.75, 0.71, 0.13, 0.07, 0.06, 0.04, 0.0, 0.0.

Table 2 :

 2 Correlation coefficients computed with (49) (the higher the better) between unsupervised and supervised metrics.

	3 :	3DFA1/GUM-EM 3DFA2/GUM-EM 3DFA1/GStudent-EM 3DFA2/GStudent-EM Mean
	3DFA1/GUM-EM	0.89	0.65	0.93	0.80	0.82
	3DFA2/GUM-EM	0.93	0.88	0.95	0.93	0.92
	3DFA3/GUM-EM	0.98	0.96	0.98	0.98	0.98
	3DFA4/GUM-EM	0.11	0.06	0.16	0.12	0.11
	3DFA1/GStudent-EM	0.80	0.57	0.88	0.74	0.75
	3DFA2/GStudent-EM	0.84	0.76	0.90	0.88	0.84
	3DFA3/GStudent-EM	0.93	0.85	0.96	0.94	0.92
	3DFA4/GStudent-EM	0.18	0.12	0.23	0.18	0.18
	GT/GUM-EM	0.73	0.54	0.82	0.71	0.70
	GT/GStudent-EM	0.67	0.48	0.78	0.66	0.65
	Alignment method		Statistical face models trained with D 1 and D 2 datasets:
	using dataset D 3 :	3DFA1/GUM-EM 3DFA2/GUM-EM 3DFA1/GStudent-EM 3DFA2/GStudent-EM Mean
	3DFA1/GUM-EM	0.30	0.28	0.28	0.29	0.29
	3DFA2/GUM-EM	0.33	0.37	0.30	0.33	0.33
	3DFA3/GUM-EM	0.16	0.17	0.16	0.15	0.16
	3DFA4/GUM-EM	0.05	0.04	0.05	0.04	0.04
	3DFA1/GStudent-EM	0.25	0.26	0.26	0.26	0.26
	3DFA2/GStudent-EM	0.25	0.28	0.23	0.26	0.25
	3DFA3/GStudent-EM	0.13	0.13	0.13	0.12	0.13
	3DFA4/GStudent-EM	0.05	0.04	0.07	0.03	0.13

Note that the translation vector t is evaluated outside the EM procedure.

https://docs.scipy.org/doc/scipy/reference/optimize. html

Acknowledgments. This work has been funded by the EU H2020 project #871245 SPRING and by the Multidisciplinary Institute in Artificial Intelligence (MIAI) # ANR-19-P3IA-0003.

Appendix A. Closed-Form Solution Using Unit Quaternions Consider (17) with Σ = σI. We immediately obtain the following formulas for the model parameters:

The formula for the posteriors becomes:

It is well known that a rotation matrix can be parameterized by a unit quaternion [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF]. Let R be parameterized by its axis and angle of rotation, n = (n 1 n 2 n 3 ) , n = 1 and φ. The unit quaternion parameterizing the rotation is:

q 0 q 1 q 2 q 3 ) ∈ R 4 by abuse of notation, and qq = 1. A vector a ∈ R 3 can be represented as a purely imaginary quaternion, namely a = (0 a 1 a 2 a 3 ) ∈ R 4 .

The action of a rotation onto a can be written as q * a * q, where the symbol * corresponds to the quaternion product and q is the conjugate of q, namely q = q 0iq 1jq 2kq 3 . Making use of the properties q 1 * q 2 2 = q 1 2 q 2 2 and q * q = q 2 = 1, the squared Euclidean norm in (A.2) can be successively written as:

= q M n q, (A.6) with:

Consequently, the right-hand side of (A.2) writes N n=1 q α n M n q = q N n=1 α n M n q = q Mq, where α n ≥ 0 and M n ∈ R 4×4 is semi-definite positive symmetric, i.e. (A.7), hence so is M. By constraining the minimizer to be a unit quaternion, we obtain the following minimization problem: min q Q(q) = min q q Mq + λ(1 -q q) .

(A.10)

From dQ/dq = 0 we obtain Mq = λq and by substitution in (A.10) we obtain Q(q ) = λ. Therefore, the minimization problem (A.10) is equivalent to estimating the smallest eigenvalueeigenvector pair (λ , q ) of M.