Équipes-Projets PERCEPTION et MISTIS et eMOTION

Keywords: Perception, cognition, interaction Point registration, feature matching, articulated object tracking, hand tracking, object pose, robust statistics, outlier detection, expectation maximization, EM, ICP, Gaussian mixture models, convex optimization, SDP relaxation

This paper addresses the issue of matching rigid and articulated shapes through probabilistic point registration. The problem is recast into a missing data framework where unknown correspondences are handled via mixture models. Adopting a maximum likelihood principle, we introduce an innovative EM-like algorithm, namely the Expectation Conditional Maximization for Point Registration (ECMPR) algorithm. The algorithm allows the use of general covariance matrices for the mixture model components and improves over the isotropic covariance case. We analyse in detail the associated consequences in terms of estimation of the registration parameters, and we propose an optimal method for estimating the rotational and translational parameters based on semi-definite positive relaxation. We extend rigid registration to articulated registration. Robustness is ensured by detecting and rejecting outliers through the addition of a uniform component to the Gaussian mixture model at hand. We provide an in-depth analysis of our method and we compare it both theoretically and experimentally with other robust methods for point registration.

Mots-clés : Appariement de points, suivi d'objets articulés, suivi de la main, pose, statistique robuste, rejet de outliers, EM, ICP, relaxation SDP.

1 Introduction, related work, and contributions

In image analysis and computer vision there is a long tradition of algorithms for finding an optimal alignment between two sets of points. This is referred to as the point registration (PR) problem, which is twofold: (i) Find point-to-point correspondences and (ii) estimate the transformation allowing the alignment of the two sets. Existing PR methods can be roughly divided into three categories:

The Iterative Closest Point (ICP) algorithm [START_REF] Besl | A method for registration of 3-D shapes[END_REF][START_REF] Zhang | Iterative point matching for registration of free-form curves and surfaces[END_REF] and its numerous extensions [START_REF] Rusinkiewicz | Efficient variants of the ICP algorithm[END_REF][START_REF] Fitzgibbon | Robust registration of 2D and 3D point sets[END_REF][START_REF] Chetverikov | Robust Euclidean alignment of 3D point sets: the trimmed iterative closest point algorithm[END_REF][START_REF] Sharp | Maximum-likelihood registration of range images with missing data[END_REF][START_REF] Demirdjian | Combining geometric-and view-based approaches for articulated pose estimation[END_REF][START_REF] Munderman | Accurately measuring human movement using articulated ICP with soft-joint constraints and a repository of articulated models[END_REF], soft assignment methods [START_REF] Rangarajan | The softassign procrustes matching algorithm[END_REF][START_REF] Chui | A new point matching algorithm for non-rigid registration[END_REF][START_REF] Liu | Automatic 3D free form shape matching using the graduated assignment algorithm[END_REF][START_REF]A mean field annealing approach to accurate free form shape matching[END_REF], and probabilistic methods [START_REF] Wells | Statistical approaches to feature-based object recognition[END_REF][START_REF] Chui | A feature registration framework using mixture models[END_REF][START_REF] Granger | Multi-scale EM-ICP: A fast and robust approach for surface registration[END_REF][START_REF] Myronenko | Non-rigid point set registration: Coherent point drift[END_REF][START_REF] Sofka | Simultaneous covariance driven correspondence (CDC) and transformation estimation in the expectation maximization framework[END_REF][START_REF] Jian | A robust algorithm for point set registration using mixture of Gaussians[END_REF] to cite just a few.

ICP alternates between binary point-to-point assignments and optimal estimation of the transformation parameters. Efficient versions of ICP use sampling processes, either deterministic or based on heuristics [START_REF] Rusinkiewicz | Efficient variants of the ICP algorithm[END_REF]. Sampling strategies can be cast into more elaborate outlier rejection methods such as [START_REF] Fitzgibbon | Robust registration of 2D and 3D point sets[END_REF] which applies a robust loss function to the Euclidean distance, thus yielding a non-linear version of ICP called LM-ICP. Another standard robust method is to select trimmed subsets of points through repeated random sampling, such as the TriICP algorithm proposed in [START_REF] Chetverikov | Robust Euclidean alignment of 3D point sets: the trimmed iterative closest point algorithm[END_REF]. In [START_REF] Sharp | Maximum-likelihood registration of range images with missing data[END_REF] a maximum-likelihood non-linear optimizer is bootstrapped by combining ICP with a RANSAC-like trimming method [START_REF] Meer | Robust techniques for computer vision[END_REF]. Although ICP is attractive for its efficiency, it can be easily trapped in local minima due to the strict selection of the best point-to-point assignments. This makes ICP to be particularly sensitive both to initialization and to the choice of a threshold needed to accept or to reject a match.

The closest-point strategy of ICP can be replaced by soft assignments within a continuous optimization framework [START_REF] Rangarajan | The softassign procrustes matching algorithm[END_REF][START_REF] Chui | A new point matching algorithm for non-rigid registration[END_REF]. Let m ji be the positive entries of the assignment matrix M, subject to the constraints j m ji = 1, i m ji = 1. When there is an equal number of points in the two sets, M is a doubly stochastic matrix. This introduces nonconvex constraints. In practice, the PR problem is solved using Lagrange parameters and a barrier function within a constrained optimization approach [START_REF] Rangarajan | The softassign procrustes matching algorithm[END_REF]. The RPM algorithm [START_REF] Chui | A new point matching algorithm for non-rigid registration[END_REF] extends [START_REF] Rangarajan | The softassign procrustes matching algorithm[END_REF] to deal with outliers. This is done by adding one column and one row to matrix M, say M. Several data points are allowed to be assigned to this extra column and, symmetrically, several model points may be assigned to this extra row. Therefore, the resulting algorithm must provide optimal entries for M and satisfy the constraints on M, thus providing one-to-one assignments for inliers, and many-to-one assignments for outliers, i.e., several entries are allowed to be equal to 1 in both the extra row and the extra column. As a consequence, M is not doubly stochastic anymore and hence the convergence properties as described in [START_REF] Sinkhorn | A relationship between arbitrary positive matrices and doubly stochastic matrices[END_REF] are not guaranteed in the presence of outliers.

Probabilistic point registration uses, in general, Gaussian mixture models (GMM). Indeed, one may reasonably assume that points from the first set (the data) are normally distributed around points belonging to the second set (the model). Therefore, the point-to-point assignment problem can be recast into that of estimating the parameters of a mixture. This can be done within the framework of maximum likelihood with missing data because one has to estimate the mixture parameters as well as the point-to-cluster assignments, i.e., the missing data. In this case the algorithm of choice is the expectation-maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood estimation from incomplete data via the EM algorithm (with discussion)[END_REF]. Formally, the latter replaces the maximization of the INRIA observed-data log-likelihood with the maximization of the expected completedata log-likelihood conditioned by the observations. As it will be explained in detail in this paper, there are intrinsic difficulties when one wants to cast the PR problem in the EM framework. The main topic and contribution of this paper is to propose an elegant and efficient way to do that.

In the recent past, several interesting EM-like implementations for point registration have been proposed [START_REF] Wells | Statistical approaches to feature-based object recognition[END_REF][START_REF] Chui | A feature registration framework using mixture models[END_REF][START_REF] Granger | Multi-scale EM-ICP: A fast and robust approach for surface registration[END_REF][START_REF] Myronenko | Non-rigid point set registration: Coherent point drift[END_REF]. In [START_REF] Wells | Statistical approaches to feature-based object recognition[END_REF] the posterior marginal pose estimation (PMPE) method estimates the marginalized joint posterior of alignment and correspondence over all possible correspondences. This formulation does not lead to the standard M-step of EM and, in particular, it does not allow the estimation of the covariances of the Gaussian mixture components. The complete-data posterior energy function is used in [START_REF] Chui | A feature registration framework using mixture models[END_REF]. This leads to an E-step which updates a set of continuous assignment variables which are similar but not identical to the standard posterior probabilities of assigning points to clusters [START_REF] Fraley | Model-based clustering, discriminant analysis, and density estimation[END_REF]. It also leads to an M-step which involves optimization of a nonlinear energy function which is approximated for simplification. The algorithms proposed in [START_REF] Wells | Statistical approaches to feature-based object recognition[END_REF][START_REF] Chui | A feature registration framework using mixture models[END_REF] do not lead to the true maximum-likelihood (ML) solution.

In [START_REF] Chui | A feature registration framework using mixture models[END_REF] as well as in [START_REF] Granger | Multi-scale EM-ICP: A fast and robust approach for surface registration[END_REF] and [START_REF] Myronenko | Non-rigid point set registration: Coherent point drift[END_REF] a simplified GMM is used, namely a mixture with a spherical (isotropic) covariance common to all the components. This has two important consequences. First, it significantly simplifies the estimation of the alignment parameters because the Manalanobis distance is replaced by the Euclidean distance: this allows the use of closed-form solutions to find the optimal rotation matrix [START_REF] Arun | Least-squares fitting of two 3-D point sets[END_REF][START_REF] Horn | Closed-form solution of absolute orientation using orhtonormal matrices[END_REF][START_REF]Closed-form solution of absolute orientation using unit quaternions[END_REF][START_REF] Umeyama | Least-squares estimation of transformation parameters between two point patterns[END_REF], as opposed to an iterative numerical solution as proposed in [START_REF] Williams | A multiple view 3D registration algorithm with statistical error modeling[END_REF]. Second, it allows connections between GMM, EM, and deterministic annealing [START_REF] Yuille | Statistical physics, mixture of distributions, and the EM algorithm[END_REF]: The common variance is interpreted as a temperature and its value is decreased at each step of the algorithm according to an annealing schedule [START_REF] Rangarajan | The softassign procrustes matching algorithm[END_REF][START_REF] Chui | A new point matching algorithm for non-rigid registration[END_REF][START_REF] Liu | Automatic 3D free form shape matching using the graduated assignment algorithm[END_REF][START_REF]A mean field annealing approach to accurate free form shape matching[END_REF][START_REF] Chui | A feature registration framework using mixture models[END_REF][START_REF] Granger | Multi-scale EM-ICP: A fast and robust approach for surface registration[END_REF][START_REF] Myronenko | Non-rigid point set registration: Coherent point drift[END_REF]. Nevertheless, the spherical-covariance assumption inherent to annealing has a number of drawbacks: anisotropic noise in the data is not properly handled, it does not use the full Gaussian model, and it does not fully benefit from the convergence properties of EM because it anneals the variance rather than considering it as a parameter to be estimated.

Another approach is to model each one of the two point sets by two probability distributions and to measure the dissimilarity between the two distributions [START_REF] Jian | A robust algorithm for point set registration using mixture of Gaussians[END_REF][START_REF] Luo | A unified framework for alignment and correspondence[END_REF][START_REF] Tsin | A correlation-based approach to robust point set registration[END_REF][START_REF] Wang | Simultaneous nonrigid registration of multiple point sets and atlas construction[END_REF]. For example, in [START_REF] Jian | A robust algorithm for point set registration using mixture of Gaussians[END_REF], each point set is modelled by a GMM where the number of components is chosen to be equal to the number of points. In the case of rigid registration, this is equivalent to replace the quadratic loss function with a Gaussian and to minimize the sum of these Gaussians over all possible point pairs. The Gaussian acts as a robust loss function. However there are two major drawbacks: The formulation leads to a non-linear optimization problem which must be solved under the nonconvex rigidity constraints, which require proper initialization. Second, the outliers are not explicitly modeled.

This paper has the following original contributions:

We formally cast the PR problem into the framework of maximum likelihood with missing data. We derive a maximization criterion based on the expected complete-data log-likelihood. We show that, within this context, the PR problem can be solved by an instance of the the expectation conditional maximization (ECM) algorithm. It has been proven that ECM is more broadly applicable than EM while it shares its desirable convergence properties [START_REF] Meng | Maximum likelihood estimation via the ECM algorithm: a general framework[END_REF]. In ECM, each M-step is replaced by a sequence of conditional maximization steps, or CM-steps. As it will be explained and detailed in this paper, ECM is particularly well suited for point registration because the maximization over the registration parameters cannot be carried out independently of the other parameters of the model, namely the covariances. For these reasons we propose the Expectation Conditional Maximization for Point Registration algorithm (ECMPR).

The vast majority of existing rigid point registration methods use isotropic covariances for reasons that we just explained. In the more general case of anisotropic covariances, we show that the optimization problem associated with rigid alignment cannot be solved in closed-form. The iterative numerical solution proposed in [START_REF] Williams | A multiple view 3D registration algorithm with statistical error modeling[END_REF] estimates the motion parameters without estimating the covariances. We propose and devise a novel solution to this problem which consists in transforming the nonconvex problem into a convex one using semi-definite positive (SDP) relaxation [START_REF] Lemarechal | SDP relaxations in combinatorial optimization from a Lagrangian viewpoint[END_REF]. Hence, rigid alignment in the presence of anisotropic covariance matrices is amenable to a tractable optimization problem.

We extend the rigid alignment solution just mentioned to articulated alignment. Based on the fact that the kinematic motion of an articulated object can be written as a chain of constrained rigid motions, we devise an incremental solution which iteratively applies the rigid-alignment solution just mentioned to the rigid parts of the kinematic chain. There are a few methods for aligning articulated objects via point registration. In [START_REF] Demirdjian | Combining geometric-and view-based approaches for articulated pose estimation[END_REF] ICP is first applied independently to each rigid part of the articulated object and next, the articulated constraints are enforced. The articulated ICP method of [START_REF] Munderman | Accurately measuring human movement using articulated ICP with soft-joint constraints and a repository of articulated models[END_REF] alternates between associating points from the two sets and estimating the articulated pose. The latter is done by minimizing a non-linear least-square error function which ensures that the rigid body parts are in an optimal pose while the kinematic joint constraints are only weakly satisfied. Our approach has two advantages with respect to these methods. First, rather than ICP, we use ECM which has proven convergence properties and which can handle inliers and outliers in a principled way. Second, our incremental rigid alignment formulation naturally enforces the kinematic constraints. As a consequence, these constraints hold exactly and there is no need to enforce them a posteriori. Moreover, the articulated registration method that we propose takes full advantage of the rigid point registration algorithm, which is quite different from data-pointto-object-part registration [START_REF] Kakadiaris | Model-based estimation of 3d human motion[END_REF][START_REF] Plaenkers | Articulated soft objects for multi-view shape and motion capture[END_REF][START_REF] Bregler | Twist based acquisition and tracking of animal and human kinematics[END_REF][START_REF] Knossow | Human motion tracking with a kinematic parameterization of extremal contours[END_REF][START_REF] Horaud | Human motion tracking by registering an articulated surface to 3-D points and normals[END_REF]. We note that our method is similar in spirit with [START_REF] Pellegrini | A generalization of the ICP algorithm for articulated bodies[END_REF]. However, the latter suffers from the limitations of ICP.

One important property of any point registration method is its robustness to outliers. Our method has a built-in outlier model, namely a uniform component that is added to the Gaussian mixture to account for non-Gaussian data, as suggested in [START_REF] Banfield | Model-based Gaussian and non-Gaussian clustering[END_REF]. This adds an improper uniformdensity component to the mixture. In theory it is attractive to incorporate the estimation of the uniform parameters into the EM algorithm. In practice, this requires an in-depth analysis of ML for the mixture in the presence of several parametric component models, which is an unsolved INRIA Figure 1: An illustration of the point registration method applied to the problem of aligning an articulated model of a hand to a set of 3D points. The 3D data are obtained by stereo reconstruction from an image pair (top). The hand model consists of 3D points lying on 16 hand parts (one root part, i.e., the palm, and 3 additional parts for each finger). The model contains 5 kinematic chains, each one is composed of the palm and one finger, i.e., 4 parts. Hence, the palm, or the root part, is common to all the kinematic chains. This articulated model has 27 degrees of freedom (3 translations and 3 rotations for the palm, 5 rotations for the thumb and 4 rotations for the index, middle, ring, and baby fingers). The result of the ECMPR-articulated algorithm is shown projected onto the left image (bottom-left) and as an implicit surface defined as a blending over the 16 hand parts (bottom-right).

problem [START_REF] Hennig | Breakdown points for maximum likelihood estimators of location-scale mixtures[END_REF][START_REF] Hennig | The noise component in model-based cluster analysis[END_REF]. We propose a treatment of the uniform component on the basis of considerations and properties that are specific to point registration. This modifies the expressions of the posteriors without adding any extra free parameters in the maximization step and without altering the general structure of the algorithm. Hence, the convergence proofs of ECM [START_REF] Meng | Maximum likelihood estimation via the ECM algorithm: a general framework[END_REF] and of EM [START_REF] Dempster | Maximum likelihood estimation from incomplete data via the EM algorithm (with discussion)[END_REF][START_REF] Redner | Mixture densities, maximum likelihood and the EM algorithm[END_REF][START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF] carry over in this case.

Our approach to outlier rejection differs from existing methods currently used in point registration. Non-quadratic robust loss functions are pro-posed in [START_REF] Fitzgibbon | Robust registration of 2D and 3D point sets[END_REF] and in [START_REF] Sofka | Simultaneous covariance driven correspondence (CDC) and transformation estimation in the expectation maximization framework[END_REF] but the drawback is that the optimization process can be trapped in local minima. This is not the case with our method because of the embedding of outlier rejection within EM. Other robust techniques such as RANSAC [START_REF] Sharp | Maximum-likelihood registration of range images with missing data[END_REF][START_REF] Meer | Robust techniques for computer vision[END_REF], least median of squares (LMS) [START_REF] Rousseeuw | Least median of squares regression[END_REF], or least trimmed squares (LTS) [START_REF] Chetverikov | Robust Euclidean alignment of 3D point sets: the trimmed iterative closest point algorithm[END_REF][START_REF] Rousseeuw | Positive-breakdown robust methods in computer vision[END_REF] must consider a very large number of subsets sampled from the two sets of points before a satisfactory solution can be found. Moreover, there is a risk that the two trimmed subsets which are eventually selected (a data subset and a model subset) contain outlying data that lead to a good fit. These random sampling issues are even more critical when one deals with articulated objects because several subsets of trimmed data points must be available, i.e., one trimmed subset for each rigid part.

We perform extensive experiments with both the ECMPR-rigid and ECMPRariculated point registration algorithms. We thoroughly study the behaviour of the method with respect to (i) the initial parameter values, (ii) the amount of noise added to the observed data, (iii) the presence of outliers, and (iv) the use of anisotropic covariances instead of isotropic ones. We illustrate the effectiveness of the method in the case of tracking a complex articulated object -a human hand composed of 5 kinematic chains, 16 parts, and 27 degrees of freedom, as shown in Fig. 1.

The remainder of this paper is organized as follows. In section 2 the PR problem is cast into the framework of ML and the proposed algorithm is outlined. In section 3 the expected complete-data log-likelihood is derived. In section 4 the EM algorithm for point registration, ECMPR, is formally derived. The algorithm is applied to rigid point sets (section 5) and to articulated point sets (section 6). Experimental results obtained both with simulated and real data are described in section 7.

Problem formulation 2.1 Mathematical notations

Throughout the paper, vectors will be in slanted bold style while matrices will be in bold style. We will consider two sets of 3-D points. We denote by Y = {Y j } 1≤j≤m the 3-D coordinates of a set of observed data points and by X = {X i } 1≤i≤n the 3-D coordinates of a set of model points. The model points lie on the surface of either a rigid or an articulated object. Hence, each model point may undergo either a rigid or an articulated transformation which will be denoted by µ : R 3 → R 3 . The 3-D coordinates of a transformed model point µ(X i ; Θ) are parameterized by Θ. In the case of rigid registration, the parameterization will consist of a 3×3 rotation matrix R and a 3×1 translation vector t. Hence, in this case we have:

µ(X i ; Θ) = RX i + t, Θ := {R, t}. (1)
We will refer to the parameter vector Θ as the registration parameters. Section 6 will make explicit the registration parameters in the case of articulated objects.

INRIA

A parameter overscripted by * , e.g., Θ * , denotes the optimal value of that parameter. The overscript denotes the transpose of a vector or of a matrix. X -Y 2 is the squared Euclidean distance and X -Y 2 Σ is the squared Mahalanobis distance, i.e. (X -Y ) Σ -1 (X -Y ) where Σ is a 3×3 symmetric positive definite matrix.

Point registration, maximum likelihood, and EM

In this paper we will formulate point registration as the estimation of a mixture of densities: A Gaussian mixture model (GMM) is fitted to the data set Y such that the centers of the Gaussian densities are constrained to coincide with the transformed model points µ(X i ; Θ), X i ∈ X . Therefore, each density in the mixture is characterized by a mean vector µ i and a covariance matrix Σ i .

In the standard mixture model approach both the means and the covariances are the free parameters. Here the means are parameterized by the registration parameters which enforce prior knowledge about the transformation that exists between the two sets of points. Therefore, the observed-data log-likelihood is a function of both the registration parameters and of the covariance matrices:

L(Θ, Σ 1 , . . . , Σ n |Y) = log P (Y; Θ, Σ 1 , . . . , Σ n ) (2)
The direct maximization of L over these parameters is intractable due to the presence of missing data, namely the unknown assignment of each observed data point Y j to one of the mixture's components. Let Z = {Z j }, 1 ≤ j ≤ m be these missing data which will be treated as a set of hidden random variables. Each variable Z j assigns an observed data point Y j to a model point X i , 1 ≤ j ≤ n, or to an outlier class indexed by n + 1.

Dempster, Laird & Rubin [START_REF] Dempster | Maximum likelihood estimation from incomplete data via the EM algorithm (with discussion)[END_REF] proposed to replace L with the expected complete-data log-likelihood conditioned by the observed data, where the term complete-data refers to both the observed data Y and the missing data Z, and where the expectation is taken over the missing data (or the hidden variables):

E(Θ, Σ 1 , . . . , Σ n |Y, Z) = E Z [log P (Y, Z; Θ, Σ 1 , . . . , Σ n )|Y] (3)
Expectation maximization (EM) [START_REF] Dempster | Maximum likelihood estimation from incomplete data via the EM algorithm (with discussion)[END_REF], is an iterative method for finding maximum likelihood estimates in incomplete-data problems like the one just stated. It has been proven that the EM algorithm converges to a local maximum of the expected complete-data log-likelihood (E) and that the maximization of E also maximizes the observed-data log-likelihood L [START_REF] Redner | Mixture densities, maximum likelihood and the EM algorithm[END_REF][START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF].

Proposed algorithm

As it will be explained in detail below, in the case of point registration, EM must be replaced by ECM, i.e.:

Provide initial values for the model parameters;

RR n°7114

2. E-step. Compute the posterior probabilities given the current estimates of the registration parameters (Θ q ) and of the covariance matrices Σ q = (Σ q 1 , . . . , Σ q n ):

α q ji = P (Z j = i|X j ; Θ q , Σ q )
3. CM-steps. Maximize the expectation in (3) with respect to:

(a) The registration parameters, conditioned by the current covariance matrices:

Θ q+1 = arg max Θ E Z [log P (Y, Z; Θ, Σ q )|Y] (b)
The covariance matrices conditioned by the newly estimated registration parameters:

Σ q+1 = arg max Σ E Z [log P (Y, Z; Θ q+1 , Σ)|Y]
4. Check for convergence.

Point registration and maximum likelihood with hidden variables

In order to estimate the registration parameters, one needs to find correspondences between the observed data points and the model points. These correspondences are the missing data and will be treated as hidden variables within the framework of maximum likelihood. Hence, there is a strong analogy with clustering. An observed data point Y j could be assigned either to a Gaussian cluster centered at µ(X i ; Θ), or to a uniform class defined in detail below. In section 2.2 we already briefly introduced the hidden variables Z = {Z j } 1≤j≤m which describe the assignments of the observations to clusters, or equivalently, the data-point-to-model-point correspondences. More specifically, the notation Z j = i (or Z : j → i) means that the observation Y j matches the model point X i while Z j = n + 1 means that the observation Y j is an outlier.

We also denote by p i = P (Z j = i) the prior probability that observation Y j belongs to cluster i with center µ(X i ; Θ) and by p n+1 = P (Z j = n + 1) the prior probability of observation j to be an outlier. We also denote by P (Y j |Z j = i), ∀j ∈ {1, . . . , m}, ∀i ∈ {1, . . . , n + 1} the conditional likelihood of Y j , namely the probability of Y j given its cluster assignment.

The likelihood of an observation j given its assignment to cluster i is drawn from a Gaussian distribution with mean µ(X i ; Θ) and covariance Σ i :

P (Y j |Z j = i) = N (Y j |µ(X i ; Θ), Σ i ), ∀i, 1 ≤ i ≤ n (4)
Similarly, the likelihood of an observation given its assignment to the outlier cluster is a uniform distribution over the volume V of the 3-D working space:

P (Y j |Z j = n + 1) = U(Y j |V, 0) = 1 V (5)

INRIA

Since {Z j = 1, . . . , Z j = n, Z j = n + 1} is a partition of the event space of Z j , the marginal distribution of an observation is:

P (Y j ) = n+1 i=1 p i P (Y j |Z j = i) (6)
By assuming that the observations are independent and identically distributed, the observed-data log-likelihood L, i.e., eq. ( 2) writes:

log P (Y) = m j=1 log n i=1 p i N (Y j |µ(X i ; Θ), Σ i ) + p n+1 V (7)
The observed-data log-likelihood is conditioned by the registration parameters Θ (which constrain the centers of the Gaussian clusters), by n covariance matrices Σ i , by n + 1 cluster priors p i subject to the constraint n+1 i=1 p i = 1, and by the uniform-distribution parameter V . In the next section we will discuss the choice of the priors and the parameterization of the uniform distribution in the specific context of point registration.

It will be convenient to denote the parameter set by:

Ψ = {Θ, Σ 1 , . . . , Σ n } (8)
A powerful method for finding ML solutions in the presence of hidden variables is to replace the observed-data log-likelihood with the complete-data log-likelihood and to maximize the expected complete-data log-likelihood conditioned by the observed data. The criterion to be maximized (i.e., eq. ( 3)) becomes [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]:

E(Ψ|Y, Z) = Z P (Z|Y, Ψ) log P (Y, Z; Ψ) (9)
4 Robust EM for point registration

In this section we formally derive the EM algorithm for robust point registration. We start by making explicit the posterior probabilities of the assignments conditioned by the observations when both the observed data and the model data are described by 3-D points; Using Bayes' rule we have:

α ji = P (Z j = i|Y j ) = P (Y j |Z j = i)P (Z j = i) P (Y j ) (10)
In general, EM treats the priors p i = P (Z j = i) as parameters. In the case of point registration we propose to specialize the priors as follows:

p i = p in = v V if 1 ≤ i ≤ n p out = V -nv V if i = n + 1 ( 11
)
where v = 4πr 3 /3 is the volume of a small sphere with radius r centered at a model point X i . We assume nv V . By combining eqs. ( 4), ( 5), [START_REF] Sharp | Maximum-likelihood registration of range images with missing data[END_REF], and (11) we obtain, for all i = 1 . . . n:

α ji = |Σ i | -1 2 exp -1 2 Y j -µ(X i ; Θ) 2 Σi n k=1 |Σ k | -1 2 exp -1 2 Y j -µ(X k , Θ) 2 Σ k + ∅ 3D (12)
RR n°7114

where ∅ 3D corresponds to the outlier component in the case of 3-D point registration:

∅ 3D = 1.5 √ 2πr -3 (13)
Note that there is a similar expression in the case of 2-D point registration, namely ∅ 2D = 2r -2 and that this can be generalized to any dimension. The posterior probability of an outlier is given by:

α j n+1 = 1 - n i=1 α ji (14)
Next we derive an explicit formula for E in eqs. ( 3) and [START_REF] Rangarajan | The softassign procrustes matching algorithm[END_REF]. For that purpose we expand the complete-data log-likelihood:

log P (Y, Z; Ψ) = log m j=1 P (Y j , Z j ; Ψ) = log m j=1 P (Y j |Z j ; Ψ)P (Z j ) = log m j=1 n+1 i=1 {p i P (Y j |Z j = i; Ψ)} δ iZ j
where δ izj is the Kronecker symbol defined by:

δ iZj = 1 if Z j = i 0 otherwise (15)
Therefore, eq. ( 3), i.e., E(Ψ|Y, Z) can be written as:

E Z   m j=1 n+1 i=1 δ iZj log p i + log P (Y j |Z j = i, Ψ) | Y   = m j=1 n+1 i=1 E Z [δ iZj |Y] log p i + log P (Y j |Z j = i, Ψ) (16)
where the conditional expectation of δ iZj writes

E Z [δ iZj |Y] = n+1 k=1 δ ik P (Z j = k|Y j ) = α ji (17)
By replacing the conditional probabilities with the normal and uniform distributions, i.e., (4) and ( 5), and by neglecting constant terms, i.e., terms that do not depend on Ψ, eq. ( 16) can be written as:

E(Ψ) = - 1 2 m j=1 n i=1 α ji Y j -µ(X i ; Θ) 2 Σi + log |Σ i | (18)
It was proven that the maximizer of (18) also maximizes the observed-data log-likelihood [START_REF] Demirdjian | Combining geometric-and view-based approaches for articulated pose estimation[END_REF] and that this maximization may be carried out by the EM INRIA algorithm [START_REF] Redner | Mixture densities, maximum likelihood and the EM algorithm[END_REF][START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF]. Nevertheless, there is an additional difficulty in the case of point registration. In the standard EM, the free parameters are the means and the covariances of the Gaussian mixture and the estimation of these parameters is quite straightforward. In the case of point registration, the means are constrained by the registration parameters and, moreover, the functions µ i (X i ; Θ) are complicated by the presence of the rotation matrices, as detailed in section 5.

In practice, the estimation of Θ is conditioned by the covariances. The simultaneous estimation of all the model parameters within the M-step would lead to a difficult non-linear minimization problem. Instead, we propose to minimize [START_REF] Jian | A robust algorithm for point set registration using mixture of Gaussians[END_REF] over Θ while keeping the covariance matrices constant, which leads to [START_REF] Meer | Robust techniques for computer vision[END_REF] below, and next we estimate the empirical covariances Σ i using the newly estimated registration parameters. This amounts to replace EM by ECM [START_REF] Meng | Maximum likelihood estimation via the ECM algorithm: a general framework[END_REF]. In practice we obtain two conditional minimization steps, using α (q) ij given by ( 12) and ( 14):

Θ q+1 = arg min Θ 1 2 m j=1 n i=1 α q ji Y j -µ(X i ; Θ) 2 Σ q i ( 19
)
and for all i = 1 . . . n,

Σ q+1 i = m j=1 α q ji (Y j -µ(X i ; Θ q+1 ))(Y j -µ(X i ; Θ q+1 )) m j=1 α q ji ( 20
)
It is well known (e.g. [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF][START_REF] Ingrassia | Constrained monotone EM algorithms for finite mixture of multivariate Gaussians[END_REF]) that when the mean µ i of one of the Gaussian components collapses onto a specific data point while the other data points are "infinitely" away from µ i , the entries of the corresponding covariance matrix Σ i tend to zero. Since [START_REF] Hathaway | A constrained formulation of maximum-likelihood estimation for normal mixture distributions[END_REF], the phenomenon has been well studied for Gaussian mixtures. Under suitable conditions, constrained global maximum likelihood formulations have been proposed, which present no singularities and a smaller number of spurious maxima (see [START_REF] Ingrassia | Constrained monotone EM algorithms for finite mixture of multivariate Gaussians[END_REF] and the references therein). However, in practice these studies do not always lead to efficient EM implementations. Thus, in order to avoid such degeneracies, the covariance is artificically fattened as follows. Let QDQ be the eigendecomposition of Σ i and let's replace the diagonal matrix D with D + εI. We obtain

Σ ε i = Q(D + εI)Q = Σ i + εI.
Hence, adding εI, where ε is a small positive number slightly fattens the covariance matrix without affecting its characteristics (eccentricity and orientation of the associated ellipsoid). A more theoretical analysis and other similar transformations of problematic covariance matrices are proposed in [START_REF] Ingrassia | Constrained monotone EM algorithms for finite mixture of multivariate Gaussians[END_REF] but the straightforward choice above provided satisfying results.

Alternatively, one may model all the components of the mixture with a common covariance matrix:

Σ q+1 = m j=1 n i=1 α q ji (Y j -µ(X i ; Θ q+1 ))(Y j -µ(X i ; Θ q+1 )) m j=1 n i=1 α q ji ( 21
)
When the number of data points is small, it is preferable to use ( 21), e.g., Fig. 2.

RR n°7114

Notice that ( 19) can be further simplified by introducing the virtual observation W i and its weight λ i that are assigned to a model point X i :

W i = 1 λ i m j=1 α ji Y j ( 22
)
λ i = m j=1 α ji (23)
By expanding [START_REF] Meer | Robust techniques for computer vision[END_REF], substituting the corresponding terms with ( 22) and ( 23), and by neglecting constant terms the minimizer yields a simpler expression:

Θ q+1 = arg min Θ 1 2 n i=1 λ q i W q i -µ(X i ; Θ) 2 Σ q i ( 24
)
It is worth noticing that the many-to-one assignment model developed here has a one-to-one (data-point-to-model-point) structure: The virtual observation W i (corresponding to a normalized sum over all the observations weigthed by their posteriors {Y j , α ji }, 1 ≤ j ≤ m) is assigned to the model point X i . Eq. ( 24) will facilitate the development of an optimization method for the rigid and articulated point registration problems as outlined in the next sections. Moreover, the minimization of ( 24) is computationally more efficient than the minimization of eq. ( 19) because it involves fewer terms. In this example we modeled all the components of the mixture model with a common covariance matrix (shown with ellipses), as in [START_REF] Dempster | Maximum likelihood estimation from incomplete data via the EM algorithm (with discussion)[END_REF]. The lines correspond to current data-to-model assignments. The algorithm stopped at the 35-th iteration. There are 12 data-to-model assignments and 7 data-to-outlier assignments corresponding to the ground truth. This example corresponds to the second row in Table 1. See Section 7 for more details.

INRIA 5 Rigid point registration

In this section we assume that the model points lie on a rigid object. Therefore:

Θ := (R, t) . (25)
Eq. ( 1) holds in this case and ( 24) becomes:

Θ * = arg min R,t 1 2 n i=1 λ i W i -RX i -t 2 Σi ( 26
)
Minimization with respect to the translation parameters is easily obtained by taking the derivatives of ( 26) with respect to the 3-D vector t and setting these derivatives to zero. We obtain:

t * = n i=1 λ i Σ -1 i -1 n i=1 λ i Σ -1 i (W i -RX i ) (27)
By substituting this expression in [START_REF] Umeyama | Least-squares estimation of transformation parameters between two point patterns[END_REF], we obtain:

R * = arg min R 1 2 n i=1 λ i X i R Σ -1 i RX i + 2X i R Σ -1 i t * -2X i R Σ -1 i W i -2t * Σ -1 i W i + t * Σ -1 i t * (28)
The minimization of (28) must be carried out in the presence of the orthonormality constraints associated with the rotation matrix, i.e., RR = I and |R| = +1.

Isotropic covariance model

Eq. ( 28) significantly simplifies when isotropic covariance matrices are being used, namely Σ i = σ 2 i I 3 . In this case, the criterion above has a much simpler form because the Mahalanobis distance reduces to the Euclidean distance. We obtain:

t * = n i=1 λ i σ -2 i (W i -RX i ) n i=1 λ i σ -2 i ( 29
)
and:

R * = arg min R 1 2 n i=1 λ i σ -2 i 2X i R t * -2X i R W i -2t * W i + t * t (30)
The vast majority of existing point registration methods use an isotropic covariance. The minimizer of (30) can be estimated in closed-form using one of the methods proposed in [START_REF] Arun | Least-squares fitting of two 3-D point sets[END_REF][START_REF] Horn | Closed-form solution of absolute orientation using orhtonormal matrices[END_REF][START_REF] Umeyama | Least-squares estimation of transformation parameters between two point patterns[END_REF].

Anisotropic covariance model

In this section we provide a solution for [START_REF] Yuille | Statistical physics, mixture of distributions, and the EM algorithm[END_REF] in the general case i.e. when the covariances are anisotropic. Our formulation relies on transforming (28) into a constrained quadratic optimization problem and on using semi-definite positive (SDP) relaxation to solve it, as detailed below. We denote by r the 9 × 2 vector containing the entries of the 3×3 matrix R, namely r := vec(R). We also denote by ρ the following rank-one positive symmetric matrix:

ρ := rr (31)
By developing and regrouping terms, (28) can be written as the following quadratic minimization criterion subject to orthogonality constraints:

r * = arg min 1 2 r Ar + 2b r r ∆ kl r = δ kl , k = 1, 2, 3; l = 1, 2, 3. ( 32
)
The entries of the 9×9 real symmetric matrix A and that of the 9×1 vector b are easily obtained by identification with the corresponding terms in [START_REF] Yuille | Statistical physics, mixture of distributions, and the EM algorithm[END_REF]; The entries of A and of b are derived in the Appendix. The entries of the six 9×9 matrices ∆ kl are easily obtained from the constraint RR = I.

As already outlined, one fundamental tool for solving such a constrained quadratic optimization problem is SDP relaxation [START_REF] Lemarechal | SDP relaxations in combinatorial optimization from a Lagrangian viewpoint[END_REF][START_REF] Lemarechal | Semidefinite relaxations and Lagrangian duality with application to combinatorial optimization[END_REF]. Indeed, a quadratic form such as r Ar can equivalently be written as the matrix dot-product A, rr 1 . Using the notation (31) one can rewrite [START_REF] Meng | Maximum likelihood estimation via the ECM algorithm: a general framework[END_REF]:

     (ρ * , r * ) = arg min (ρ,r) 1 2 A, ρ + 2b r ∆ kl , ρ = δ kl , k = 1, 2, 3; l = 1, 2, 3. ρ = rr (33)
In [START_REF] Lemarechal | SDP relaxations in combinatorial optimization from a Lagrangian viewpoint[END_REF] everything is linear except the last constraint which is nonconvex. As already noticed, matrix rr is a rank-one positive symmetric matrix. Relaxing the positivity constraint to semi-definite positivity amounts to taking the convex hull of the rank-one positive symmetric matrices. Within this context, [START_REF] Lemarechal | SDP relaxations in combinatorial optimization from a Lagrangian viewpoint[END_REF] relaxes to:

     (ρ * , r * ) = arg min (ρ,r) 1 2 A, ρ + 2b r ∆ kl , ρ = δ kl , k = 1, 2, 3; l = 1, 2, 3. ρ rr (34)
To summarize, rigid point registration with anisotropic covariances, i.e. [START_REF] Yuille | Statistical physics, mixture of distributions, and the EM algorithm[END_REF], can be formulated as the convex optimization problem [START_REF] Kakadiaris | Model-based estimation of 3d human motion[END_REF]. It is well known that this generally provides a very good initial solution to a standard non-linear optimizer such as the one proposed in [START_REF] Williams | A multiple view 3D registration algorithm with statistical error modeling[END_REF]. Finally, this yields the following algorithm illustrated in Fig. 2:

1 The dot-product of two n × n matrices A = (A ij ) i,j and B = (B ij ) i,j is defined as

A, B := n P i=1 n P j=1 A ij B ij .

INRIA

The ECMPR-rigid algorithm:

1. Initialization. Set R q = I, t q = 0. Choose the initial covariance matrices Σ q i , i = 1 . . . n.

2. E-step. Evaluate the posteriors α q ji from ( 12) and ( 14), W q i from ( 22), and λ q i from ( 23), using the current parameters R q , t q , and Σ q i .

3. CM-steps.

(a) Use SDP relaxation to estimate the new rotation matrix R q+1 by minimization of ( 28) with the current posteriors α q ji and the current covariances Σ q i ; (b) Estimate the new translation vector t q+1 from (27) using the new rotation R q+1 and the current posteriors and covariance matrices;

(c) Estimate the new covariances from [START_REF] Sinkhorn | A relationship between arbitrary positive matrices and doubly stochastic matrices[END_REF] or from [START_REF] Dempster | Maximum likelihood estimation from incomplete data via the EM algorithm (with discussion)[END_REF] with the current posteriors and the new rotation matrix and translation vector.

4.

Convergence. Compare the new and current rotations. If R q+1 -R q 2 < ε then go to the Classification step. Else, set the current parameter values to their new values and return to the E-step.

5.

Classification. Assign each observation to a model point (inlier) or to the uniform class (outlier) based on the maximum a posteriori (MAP) principle:

z j = arg max i α q ji .
6 Articulated point registration

The kinematic model

In this section we will develop a solution for the articulated point registration problem. We will consider the case of an open kinematic chain. Such a chain is generally composed of rigid parts. Two adjacent parts are mechanically linked. Each link has one, two or three rotational degrees of freedom. i.e., spherical motions. In addition we assume that the root part of such an open chain may undergo a free motion with six degrees of freedom. Consequently the articulated object motions considered here are combinations of free and constrained motions. This is more general than traditional open or closed kinematic chains considered in standard robotics [START_REF] Mccarthy | Introduction to Theoretical Kinematics[END_REF].

More precisely, any rigid part p, 1 ≤ p ≤ P , moves with respect to the root part p = 0 through a chain of constrained motions. The root part itself under- goes a free motion with up to six degrees of freedom, three rotations and three translations. We assume that a partition of the set of model points is provided, X = {X 0 , . . . X p , . . . X P }; Each subset of model points X p = {X (p) i }, 1 ≤ i ≤ n p is attached to the p th rigid part of the articulated object. It is worthwhile to point out that a partitioning of the set of observations is not required in advance and is merely an output of our method. The model point X (p) i belonging to part p is transformed with:

RR n°7114 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k)
µ(X (p) i ; Θ) = R(Θ)X (p) i + t(Θ), Θ := {Θ 0 , . . . Θ p } (35)
The main difference between rigid and articulated motion is that in the former case, (i.e., eq. ( 1)) the rotation matrix and translation vector are the free parameters while in the latter case, (i.e., eq. ( 35)) the motion of any part is constrained by both the kinematic parameters Θ 1 , . . . Θ p and by the motion of the root part Θ 0 . For convenience we will adopt the homogeneous representation of the Euclidean group of 3-D rigid displacements. Hence the rotation matrix and the translation vector can be embedded into a 4×4 displacement matrix T p (Θ). The latter may well be written as a chain of homogeneous transformations:

T p (Θ) = Q 0 (Θ 0 )Q 1 (Θ 1 ) . . . Q p (Θ p ) (36)
Q 0 describes the free motion of the root part parameterized by Θ 0 = {vec(R 0 ), t 0 }.

INRIA

Each transformation Q p , 1 ≤ p ≤ P has two components: A fixed component that describes a change of coordinates, and a constrained motion component parameterized by one, two or three angles [START_REF] Knossow | Human motion tracking with a kinematic parameterization of extremal contours[END_REF].

Therefore, the estimation of the parameter vector Θ amounts to solving a difficult inverse kinematic problem, namely a set of non-linear equations that are generally solved using iterative optimization methods requiring proper initialization.

Rather than estimating the problem parameters simultaneously, in this section we devise a closed-form solution which is based on the formulation developed in section 5. We propose the ECMPR-articulated algorithm which is built on top of the ECMPR-rigid algorithm and which solves for the free and kinematic parameters incrementally by considering a single rigid part at each iteration. This contrasts with methods that attempt to estimate all the kinematic parameters simultaneously from data-point-to-object-part associations, as done in previous approaches [START_REF] Kakadiaris | Model-based estimation of 3d human motion[END_REF][START_REF] Plaenkers | Articulated soft objects for multi-view shape and motion capture[END_REF][START_REF] Bregler | Twist based acquisition and tracking of animal and human kinematics[END_REF][START_REF] Knossow | Human motion tracking with a kinematic parameterization of extremal contours[END_REF][START_REF] Horaud | Human motion tracking by registering an articulated surface to 3-D points and normals[END_REF]. The motion of the root part of the articulated object is parameterized by a rotation and a translation, while the motion of each one of the other parts is parameterized by a rotation, hence Θ 0 = {vec(R 0 ), t 0 } and Θ p = {vec(R p )}, for all p = 1 . . . P , and:

Q 0 = R 0 t 0 0 1 , Q p = R p 0 0 1 . (37)
Moreover, eq. ( 36) can be written as T 0 = Q 0 , T 1 = T 0 Q 1 , or more generally ∀p, 1 ≤ p ≤ P :

T p = T p-1 Q p ( 38
)
which can be expanded as:

T p = R 0,p-1 R p t 0,p-1 0 1 ( 39
)
where the rotation matrix R 0,p-1 and the translation vector t 0,p-1 are associated with T p-1 describing the articulated pose of body part p -1.

The pose of an articulated shape

As already mentioned, there are n p model points X (p) i associated with the p th body part. Using the set of available observations together with current estimates of their posterior probabilities one can easily compute the set of n p virtual observations W i and their weights λ i . Therefore, the criterion [START_REF] Umeyama | Least-squares estimation of transformation parameters between two point patterns[END_REF] allows rigid registration of the root part as well as registration of the p th body part conditioned by the articulated pose of the (p -1) th body part:

Θ * 0 = arg min R0,t0 1 2 n0 i=1 λ i W i -R 0 X (0) i -t 0 2 Σi (40) R * p = arg min Rp 1 2 np i=1 λ i W i -R 0,p-1 R p X (p) i -t 0,p-1 2 Σi ( 41
)

RR n°7114

By introducing the following substitutions:

U p = R 0,p-1 R p R 0,p-1 (42)
V (p) i = R 0,p-1 X (p) i ( 43
)
the minimization of ( 41) becomes:

U * p = arg min Up 1 2 np i=1 λ i W i -U p V (p) i -t 0,p-1 2 Σi ( 44
)
Therefore, if the transformation T p-1 is known, the parameters of the transformation Q p can be obtained by minimization of ( 44) and from ( 42) and ( 43), which is strictly equivalent to the minimization of [START_REF] Umeyama | Least-squares estimation of transformation parameters between two point patterns[END_REF]. To summarize, we obtain the following algorithm illustrated in Fig. 3:

The ECMPR-articulated algorithm:

1. Rigid registration of the root part. Set the current set of data points Y (0) to the whole data set. Apply the ECMPR-rigid algorithm to the data set Y (0) and to the set of model points associated with the root part (X 0 ) to estimate the pose of the root part. Compute T 0 using [START_REF] Knossow | Human motion tracking with a kinematic parameterization of extremal contours[END_REF]. Classify the data points into inliers and outliers. Remove the inliers from Y (0) to generate a new data set Y (1) .

2.

For each p = 1 . . . P : Rigid registration of the p th part. Apply the ECMPR-rigid algorithm to the current set of data points Y (q) and to the set X p . Estimate R p from ( 42) and [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF]. Compute Q p and then T p using ( 37) and [START_REF] Horaud | Human motion tracking by registering an articulated surface to 3-D points and normals[END_REF].

Classify the data points into inliers and outliers. Remove the inliers from Y (q) to generate Y (q+1) .

Experimental results

We carried out a large number of experiments with both algorithms. ECMPRrigid was applied to simulated data to assess the performance of the method with respect to (i) the initialization of the method's parameters, (ii) the amount of Gaussian noise added to the data, and (iii) the percentage of outliers drawn from a uniform distribution. ECMPR-articulated was applied to a simulated data set to illustrate the method, Fig. 3, and to the problem of tracking an articulated model of a hand using both simulated and real data.

Experiments with ECMPR-rigid

We carried out several experiments with ECMPR-rigid and with the Trimmed Iterative Closest Point algorithm (TriICP) [START_REF] Chetverikov | Robust Euclidean alignment of 3D point sets: the trimmed iterative closest point algorithm[END_REF], which is a robust implementation of ICP using random sampling. These experiments are summarized in Table 1 and on Fig. 2 and Fig. 4. In all these experiments we considered 15 model In the examples shown in Table 1 and on Fig. 2 the inliers are rotated by 25 0 and then translated using a randomized vector. The first example (first row in Table 1) is noise free. We simulated anisotropic Gaussian noise that was added to the inliers in the second and third examples. This noise is centered at each inlier location and is drawn from two one-dimensional Gaussian probability distributions with two different variances along each dimension. The variances were allowed to vary between 10% and 100% of box bounding the set of observations. In all the reported experiments, ECMPR-rigid was initialized with a null rotation angle (the identity matrix), a null translation vector, and with large variances. We used the same data with TriICP. Unlike our method, ICP methods require proper initialization, in particular in the presence of outliers.

TriICP embeds multiple initializations using a random sampling strategy, which explains the large number of iterations of this method [START_REF] Chetverikov | Robust Euclidean alignment of 3D point sets: the trimmed iterative closest point algorithm[END_REF].

Additionally, we performed a large number of trials with ECMPR-rigid in the anisotropic covariance case (second row in Table 1). The inliers are rotated with an angle that varies between 0 0 and 180 0 . For each angle we performed 1,000 trials. Fig. 4 shows the percentage of correct matches (a), the relative error in rotation (b), and the relative error in translation (c) as a function of the rotation angle between the set of data and model points. The plotted curves correspond to the mean values and to the variances computed over 1,000 trials for each rotation.

ECMPR-rigid behaves very well in the presence of both high-amplitude anisotropic Gaussian noise and outliers. The anisotropic covariance model advocated in this paper yields better results than the isotropic model both in terms of parameter estimation and number of correct assignments. The errors in rotation and translation are consistent with the level of noise added to the inliers; Overall the performance of ECMPR-rigid is very robust in the presence of outliers. This is a crucial feature of the algorithm that directly conditions the robustness of ECMPR-articulated, since the former resides in the inner loop of the latter.

To farther assess the algorithms' performance, we computed the percentage of correct matches (see Table 1), namely the number of observations that were correctly classified over the total number of observations. In case of ECMPRrigid, this classification is based on the maximum a posteriori (MAP) principle: each observation j is assigned to the cluster i (either a Gaussian cluster for a model point or a uniform cluster for an outlier) for which the posterior probability α ji is maximum. This implies that each data point, which is not an outlier, is assigned to one model point but there may be several data points assigned to the same model point. ICP algorithms use a different assignment strategy, namely they retain the closest data point for each model point and they apply a threshold to this point-to-point distance to decide whether the assignment should be validated or not. For these reasons, the counting of matches has a different meaning with ECMPR and with ICP. For example, in the case of an anisotropic covariance model (Table 1 second row and Fig. 2), ECMPR assigned 3 outliers to 3 model points while 3 inliers were incorrectly assigned. In the case of an isotropic covariance model (Table 1, third row), 4 outliers were assigned to 4 model points while 8 inliers were incorrectly assigned. In the presence of both anisotropic noise and outliers, TriICP rejected 18 data points, namely 10 outliers and 8 inliers. Comparing correct matches then may not be straightforward. A more meaningful comparison can be made by looking at the transformation estimation. It appears that ECMPR has superior performance with smaller rotation and translation errors.

As we already mentioned and as observed by others, the initialization of TriICP (and more generally of ICP algorithms) is crucial to obtain a good match. Starting from any initial guess, ICP converges very fast (4 to 5 iterations on average). However, ICP is easily trapped in a local minimum. To overcome this problem, TriICP combines ICP with a random sampling method: The space of rotational parameters is uniformly discretized and an initial solution is randomly drawn from this space. In this example, ECMPR found 95 inliers while ICP found 177 inliers. Both algorithms estimated the correct rotation. ICP failed to estimate the correct translation (see Table 2 for a quantitative comparison). We also applied both ECMPR-rigid and ICP to real data obtained with a stereo camera pair as shown on Fig. 5 and Table 2: The two stereo image pairs of a walking person were grabbed at two diffierent time instances. Two sets of 3D points were reconstructed from these two image pairs, Fig. 5-(c). The first set has 223 "model" points and the second set has 249 "data" points. These 3D points belong either to the walking person or to the static background. Fig. 5-(d) shows the matches found by ECMPR-rigid and Fig. 5-(e) shows the matches found by ICP. Table 2 summarizes the results. Both algorithms were initialized with R = I and t = 0. The error in translation is computed with tt g / t g where t is the estimated translation vector and t g is the ground truth. The minimization error is computed with the square root of 1/n in nin i=1 Y i -RX it 2 where n in is the number of inliers estimated by each algorithm. ICP was run with different threshold values. In all cases (ECMPR and ICP) the rotation matrix is correctly estimated.

Experiments with ECMPR-articulated

We tested ECMPR-articulated on a hand tracking task, with both simulated and real data. The hand model used in all our experiments consists in five kinematic chains that share a common root part -the palm. Each kinematic chain is composed of four rigid parts, one part for the palm and three other parts for the phalanges composing each finger. Altogether, the kinematic hand model has 16 rigid parts and 21 rotational degrees of freedom (5 rotations for the thumb and 4 rotations for the other fingers). With the additional six degrees of freedom (three rotations and three translations) associated with the free motion of the palm, the hand has a total of 27 degrees of freedom. Each hand-part is modeled with an ellipsoid with fixed dimensions. Model points are obtained by uniformly sampling the surface of each one of these ellipsoids. This representation also allows to define an articulated implicit surface over the set of ellipsoids [START_REF] Plaenkers | Articulated soft objects for multi-view shape and motion capture[END_REF][START_REF] Horaud | Human motion tracking by registering an articulated surface to 3-D points and normals[END_REF][START_REF] Dewaele | Hand motion from 3D point trajectories and a smooth surface model[END_REF]. Here we only use this implicit surface representation for visualization purposes.

In the case of simulated data, we animated this model in order to produce realistic articulated motions and to generate sets of model points, one set for each pose of the model. In practice, all the experiments described below used 15 model points for each hand part which corresponds to a total of 240 model points namely X (p) i with 1 ≤ i ≤ 15 and 0 ≤ p ≤ 15. In order to simulate realistic observations we added Gaussian noise to the surface points. The standard deviation of the noise was 10% of the size of the bounding box of the data set. We also added outliers drawn from a uniform distribution defined over the volume occupied by the working space of the hand. In all these simulations the data sets contain 30% of outliers, i.e., there are 240 model points, 240 inliers and 72 outliers. Fig. 6 and Fig. 8 show two experiments performed with simulated hand motions. Each one of these simulated data (top rows) contains a sequence of 120 articulated poses. We applied our registration method to these sequences, we estimated the kinematic parameters, and we compared them with the ground truth. ECMPR-articulated is applied in parallel to the five kinematic chains. First, ECMPR-rigid registers the root part (the hand palm) common to all the chains. Second, ECMPR-rigid is applied to the first phalanx of the index, middle, ring, and baby fingers. Third, it is applied to the second phalanx, etc. Fig. 6 shows a sequence of simulated poses (top row) and the results obtained with our algorithm (middle and bottom rows). When starting with a large covariance, ECMPR correctly estimated the articulated poses of the simulated hand (middle row). Starting with small covariances is equivalent to consider the data points that are in the neighbourhood of the model points and to disregard data points that are farther away from the current model point positions. In this case the trajectory of the thumb has been correctly estimated but the other four fingers failed to bend (bottom row). Notice, however, that in both cases the tracker has been able to "catch up" with these finger motions and to reduce the discrepancy between the estimated trajectories and the ground truth. The simulated trajectories and the estimated trajectories of the first and These experiments yielded very good results. As expected, the percentage of outliers barely affected the registration results. These experiments confirmed the importance of using an anisotropic covariance model as well as the fact that covariance initialization is crucial. All the instances of the ECMPR-rigid algorithm (embedded in ECMPR-articulated) are initialized with large spherical covariances. While this increases the number of EM iterations, it allows the algorithm to escape from local minima.

INRIA

We then tested our method with real data consisting in several hand motions observed with a stereoscopic camera system, Fig 1 . Each data sequence that we used contains 100 image pairs gathered at 20 frames per second. We run a standard stereo algorithm to estimate 3-D points. This yielded 500 to 1000 reconstructed points at each time step. The noise associated with these stereo data is inherently anisotropic because of the inaccuracy in depth. Moreover, there are many outliers that correspond either to data points which do not lie on the hand or to stereo mismatches.

The results of applying ECMPR to these data sets are illustrated on Figs. 9 and 10. In the first example the hand performs a grasping movement. In the second example the hand rotates around an axis which is roughly parallel to the image plane. In both cases the algorithm selected, on average, 250 inliers per frame; This number roughly corresponds to the number of model points being considered (240). All the other data points were assigned to the outlier class. Notice that the number of data points vary a lot (500 to 1000 observations at each frame) and that the outlier rejection mechanism that we propose in this paper does not need to know in advance the percentage of outliers.

Note that along these motion sequences the hand flips from one side to another side while the positions and orientations of the fingers vary considerably. This means that it is often the case that almost all the model points that were currently registered, may suddenly disappear while other model points appear. This is one of the main difficulties associated with registering articulated objects. Therefore, during the tracking, the algorithm must perform some form of bootstrapping, i.e., it must establish data-point-to-model-point assignments Figure 9: The image of a hand and the result of tracking for a grasping movement. Notice the weird position of the thumb due to the fact that the kinematic model that we use for the thumb is not anatomically correct. from scratch. Re-initialization of the covariance matrix at each time step, along the lines described above, is crucial to the success of the registration/tracking algorithm.

INRIA

Conclusions

In this paper we addressed the problem of matching rigid and articulated shapes through robust point registration. The proposed approach has its roots in model-based clustering [START_REF] Fraley | Model-based clustering, discriminant analysis, and density estimation[END_REF]. More specifically, the point registration problem is cast into the framework of maximum likelihood with hidden variables [START_REF] Dempster | Maximum likelihood estimation from incomplete data via the EM algorithm (with discussion)[END_REF][START_REF] Redner | Mixture densities, maximum likelihood and the EM algorithm[END_REF][START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF]. We formally derived a variant of the EM algorithm which maximizes the expected complete-data log-likelihood. This guarantees maximization of the observeddata log-likelihood. We showed that it is convenient to replace the standard M-step by three conditional maximization steps, or CM-steps, while preserving the convergence properties of EM.

Our approach differs significantly from existing methods for point registration, namely ICP and its variants [START_REF] Besl | A method for registration of 3-D shapes[END_REF][START_REF] Zhang | Iterative point matching for registration of free-form curves and surfaces[END_REF][START_REF] Rusinkiewicz | Efficient variants of the ICP algorithm[END_REF][START_REF] Fitzgibbon | Robust registration of 2D and 3D point sets[END_REF][START_REF] Chetverikov | Robust Euclidean alignment of 3D point sets: the trimmed iterative closest point algorithm[END_REF][START_REF] Sharp | Maximum-likelihood registration of range images with missing data[END_REF][START_REF] Demirdjian | Combining geometric-and view-based approaches for articulated pose estimation[END_REF][START_REF] Munderman | Accurately measuring human movement using articulated ICP with soft-joint constraints and a repository of articulated models[END_REF], soft assignment methods [START_REF] Rangarajan | The softassign procrustes matching algorithm[END_REF][START_REF] Chui | A new point matching algorithm for non-rigid registration[END_REF][START_REF] Liu | Automatic 3D free form shape matching using the graduated assignment algorithm[END_REF][START_REF]A mean field annealing approach to accurate free form shape matching[END_REF], as well as various EM implementations [START_REF] Wells | Statistical approaches to feature-based object recognition[END_REF][START_REF] Chui | A feature registration framework using mixture models[END_REF][START_REF] Granger | Multi-scale EM-ICP: A fast and robust approach for surface registration[END_REF][START_REF] Myronenko | Non-rigid point set registration: Coherent point drift[END_REF][START_REF] Sofka | Simultaneous covariance driven correspondence (CDC) and transformation estimation in the expectation maximization framework[END_REF][START_REF] Jian | A robust algorithm for point set registration using mixture of Gaussians[END_REF]: The ECMPR-rigid and -articulated algorithms that we proposed fit a set of model points to a set of data points where each model point is the center of a Gaussian component in a mixture model. Each component in the mixture may have its own anisotropic covariance. Our method treats the data points and the model points in a non-symmetric way, which has several advantages: It allows to deal with a varying number of observations, either larger or smaller than the number of model points, it performs robust parameter estimation in the presence of data corrupted with noise and outliers, and it is based on a principled probabilistic approach.

More specifically, the method guarantees robustness via a uniform component added to the Gaussian mixture model. This built-in outlier rejection mechanism differs from existing outliers detection/rejection strategies used in conjunction with point registration, such as methods based on non-linear loss functions that can be trapped in local minima, or methods based on random sampling which are time-consuming and that can only deal with a limited number of outlying data.

In particular we put emphasis on a general model that uses anistropic covariance matrices, in which case the rotation associated with rigid alignment cannot be found in closed-form. This led us to approximate the associated non-convex optimization problem with a convex one. Namely, we showed how to transform the non-linear problem into a constrained quadratic optimization one and how to use semi-definite positive relaxation to solve it in practice.

We provided in detail the ECMPR-rigid algorithm. We showed how this algorithm can be incrementally applied to articulated registration using a novel kinematic representation that is well suited in the case of point registration.

In general, ECMPR performs better than ICP. In particular it is less sensitive to initialization and it is more robust to outliers. In the future we plan to investigate various ways of implementing our algorithm more efficiently. Promising approaches are based on modifying the standard E-step. A fast but suboptimal INRIA "winner take all" variant is Classification EM, or CEM, which consists in forcing the posterior probabilities to either 0 or 1 after each E-step [START_REF] Celeux | A classification EM algorithm for clustering and two stochastic versions[END_REF]. We plan to study CEM in the particular context of point registration and, possibly, derive a more efficient implementation of ECMPR. This may also lead to a probabilistic interpretation of ICP, and hence to a better understanding of the links existing between probabilistic and deterministic registration methods. Other efficient variants of the E-step are based on structuring the data using either block-like organizations [START_REF] Thiesson | Accelerating EM for large databases[END_REF], or KD-trees [START_REF] Verbeek | Accelerated em-based clustering of large data sets[END_REF]. We also plan to implement KD-trees in order to increase the efficiency of ECMPR.

A Expansion of A and b in eq. [START_REF] Meng | Maximum likelihood estimation via the ECM algorithm: a general framework[END_REF] By expanding [START_REF] Yuille | Statistical physics, mixture of distributions, and the EM algorithm[END_REF], substituting the optimal translation with ( 27) and rearranging terms, one obtains the following expressions for the 9×9 matrix A and the 9×1 vector b:

A = N -M KM (45) b = M p -q (46)
with:

N 9×9 = n i=1 λ i X i X i ⊗ Σ -1 i M 3×9 = n i=1 λ i X i ⊗ Σ -1 i K 3×3 = n i=1 λ i Σ -1 i -1 p 3×1 = K n i=1 λ i Σ -1 i W i q 9×1 = vec n i=1 λ i Σ -1 i W i X i
The Kronecker product between the m × n matrix/vector A and the p × q matrix/vector B is the mp × nq matrix/vector defined by:

A ⊗ B =

Figure 2 :

 2 Figure2: Illustration of the ECMPR-rigid algorithm. There are 15 model points (filled circles) and 25 data points, 15 inliers (filled squares) that correspond to model points that were rotated, translated and corrupted by anisotropic Gaussian noise, and 10 outliers (empty squares) drawn from a uniform distribution. In this example we modeled all the components of the mixture model with a common covariance matrix (shown with ellipses), as in[START_REF] Dempster | Maximum likelihood estimation from incomplete data via the EM algorithm (with discussion)[END_REF]. The lines correspond to current data-to-model assignments. The algorithm stopped at the 35-th iteration. There are 12 data-to-model assignments and 7 data-to-outlier assignments corresponding to the ground truth. This example corresponds to the second row in Table1. See Section 7 for more details.

Figure 3 :

 3 Figure 3: Illustration of the ECMPR-articulated algorithm. The articulated object consists of 4 rigid parts. (a)-(c): 1st, 4th and 25th iterations of ECMPRrigid used to register the root part, i.e., p = 0; the algorithm converged in 25 iterations. (d)-(f): The unmatched data (the outliers), are used to register the second part, p = 1, in 13 iterations. (g)-(i): The third part, p = 2, is aligned with the remaining data after 13 iterations. (j)-(k): The fourth part, p = 3, is registered with the remaining data in 4 iterations.

Figure 4 :

 4 Figure 4: Statistics obtained with ECMPR-rigid over a large number of trials. The percentage of correct matches, and relative errors in rotation and translation are shown as a function of the rotation angle between the data and the model points, in the presence of outliers. The plots correspond to the means (central curves) and to the means +/-the standard deviation (upper and lower curves) computed over 1,000 trials.

Figure 5 :

 5 Figure 5: Comparison between ECMPR-rigid and ICP applied to stereo data. (a) The first stereo image pair of a walking person. (b) The second stereo pair. The person performed a translational motion of 280 mm towards the camera and from right to left. (c) The two sets of 3D points before registration (223 model points and 249 data points). The result of matching with (d) ECMPR and with (e) ICP are shown superimposed onto the left image of the first pair.In this example, ECMPR found 95 inliers while ICP found 177 inliers. Both algorithms estimated the correct rotation. ICP failed to estimate the correct translation (see Table2for a quantitative comparison).

Figure 6 :

 6 Figure 6: Top row: the ground-truth of the simulated poses and the simulated data (inliers and outliers). Milddle row: a correct registration result. Bottom row: ECMPR failed to correctly estimate all the kinematic parameters due to an improper initialization of the covariance matrix.

Figure 7 :Figure 8 :

 78 Figure 7: (a) Simulated trajectories of two angular joints associated with the first and second phalanges of the index finger. (b) Good estimation of these angle values. (c) Bad estimation of the angle values due to improper initialisation of the covariance matrix in ECMPR. These trajectories correspond to the examples shown in Fig. 6.

Figure 10 :

 10 Figure 10: A rotational movement of the hand around an axis parallel to the image plane causes the fingers to disappear from the left-hand side of the image and appear again onto the right-hand side. These occlusions have as a result a very coarse initialization of the current pose. In spite of this problem the tracker performs quite well due to re-initialization of the covariance matrix at each time step of the tracker.

  (A) returns the mn × 1 vector:vec(A) = (A 11 . . . A mn ) RR n°7114

  Dans cet article on s'intèrese au problème d'appariement rigide et articulé utilisant une méthode probabiliste. Le problème est vu dans le cadre des modèles de mélanges avec variables cachées. Nous introduisons un algorithme du type EM, Expectation Conditional Maximization for Point Registration (ECMPR). Ce nouvel algorithme permet l'utilisation de covariances anisotropiques pour les composantes du modèle, en améliorant ainsi le cas des covariances isotropiques. Nous analysons en détail les conséqunces sur l'estimation des paramètres d'appariement et nous proposons une méthode optimale pour estimer la rotation et la translation utilisant la méthode de semi-definite positive relaxation. Nous étendonds le résultat au cas d'objets articulés. La détection et le rejet d'outliers assurent la robustesse de la méthode. Finalement, on compare notre méthode avec d'autres méthodes d'appariement robuste de points.

	Appariement Rigide et Articulé de Points avec
	l'Algorithme Expectation Conditional
	Maximization

Résumé :

Table 1 :

 1 Summary of experiments with simulated data using ECMPR-rigid and TriICP. the clusters' centers in the mixture model, as well as 25 observations: 15 inliers and 10 outliers. The inliers are generated from the model points: they are rotated, translated, and corrupted by noise. All the outliers in all the experiments are drawn from a uniform distribution spanning the bounding box of the set of observations.

	Method	Sim.	Cov.	No. of	Error in	Error in	Correct	Proc.
		noise	model	iter.	rotation	translation	matches	time
					(%)	(%)	(%)	(ms)
	ECMPR -	anisotr. 20	0.0	0.0	100	16.8
	ECMPR	anisotr. anisotr. 36	1.5	5.6	76	31.4
	(Fig. 2)
	ECMPR anisotr. isotr.	35	8.1	26.3	52	16.8
	TriICP	-	-	217	0.0	0.0	100	67.2
	TriICP	anisotr. -	215	10.3	6.5	28	63.2
	points corresponding to

Table 2 :

 2 Comparison between ECMPR-rigid and ICP applied to the stereo data of Fig.5.

	Algorithm	Number of	Number of	Translation	Minimization
		iterations	inliers	error (%)	error (mm)
	ECMPR	12	95	12.5	4.6
	ICP (Worst) 4	227	34.6	6.0
	ICP (Best)	6	177	22.4	4.9

RR n°7114

http://www.inria.fr ISSN 0249-6399


Mathieu Lefebvre
email: mathieu.lefebvre@ulg.ac.be

Ferdinand M Vieider #
email: vieider@gate.cnrs.fr

Marie Claire Villeval
email: villeval@gate.cnrs.fr

Incentive Effects on Risk Attitude in Small Probability Prospects

Keywords: C91, D81, D89 Risk attitude, Incentives, Decision, Experiment

Most studies on the role of incentives on risk attitude report data obtained from within-subject experimental investigations. This may however raise an issue of sequentiality of effects as later choices may be influenced by earlier ones. This paper reports instead between-subject results on the effect of monetary stakes on risk attitudes for small probability prospects in a laboratory experiment. Under low stakes, we find the typical risk seeking behavior for small probabilities predicted by the prospect theory. But under high stakes, we provide some evidence that risk seeking behavior is dramatically reduced. This could suggest that utility is not consistently concave over the outcome space, but rather contains a convex section for very small amounts.

INTRODUCTION

Notwithstanding a large research effort in individual decision making under risk, certain voids remain in our understanding of risk attitudes. One issue that has been hotly debated in the literature is the effect-and for that matter, the necessity-of the provision of monetary incentives when studying risk attitudes. As for many other economic decisions (Camerer & Hogarth, 1999;[START_REF] Johansson-Stenman | Does Stake Size Matter in Trust Games?[END_REF][START_REF] Kocher | Does Stake Size Matter for Cooperation and Punishment?[END_REF], the effect of incentives is a potentially contentious issue, since many of the traditional findings on risk attitudes have been obtained with hypothetical payoffs. After many years of heated debate, a consensus on these issues seems to be emerging. Indeed, financial incentives are generally thought to leave the qualitative findings obtained with hypothetical studies intact [START_REF] Battalio | Testing Between Alternative Models of Choice under Uncertainty[END_REF]. Quantitatively however, incentives seem to matter inasmuch as higher stakes increase risk aversion (Astrebo, Mata, & Santos-Pinto, 2009;[START_REF] Binswanger | Attitudes toward Risk: Experimental Measurement in Rural India[END_REF]Kachelmeier & Shahata, 1992). Also, while the size of real stakes matters, so do the nominal stakes in hypothetical choices [START_REF] Kühberger | Framing Decisions: Hypothetical and Real[END_REF][START_REF] Holt | Risk Aversion and Incentive Effects[END_REF]Laury & Holt, 2005).

Even though the studies cited above have accumulated a substantial (and generally coherent) view on the issue, some methodological doubts remain. Indeed, although the studies cited used different methodologies for the elicitation of risk attitudes and thus proved the stability of the finding, they all report data obtained from within-subject investigations of the issue. While within-subject investigations are statistically powerful and avoid potential confounds, they also pose an issue of sequentiality of effects as later choices may be influenced by earlier ones. For instance, Kachelmeier & Shehata remark how "the transparent manipulation of prize level may have acted as a cue to subjects that their responses should change" (p. 1131). Read (2005) criticizes [START_REF] Holt | Risk Aversion and Incentive Effects[END_REF] for repeatedly telling subjects that certain choices were hypothetical, and especially for emphasizing the contrast of those hypothetical choices to preceding and following real choices. More generally, the withinsubject method may suffer from experimenter demand effects, whereby subjects try to conform to the presumed expectations of the experimenter [START_REF] Minor | Experimenter-Expectancy Effect as a Function of Evaluation Apprehension[END_REF][START_REF] Sawyer | Demand Artifacts in Laboratory Experiments in Consumer Research[END_REF]. Even in the absence of a motivation to please the experimenter, variations could be introduced by issues of sensitization [START_REF] Greenwald | Within-Subject Designs: To Use or Not to Use[END_REF].

To the best of our knowledge, no systematic investigation of between-subject effects of different stakes on risk attitude exists. This is all the more surprising since the papers cited above seem all to have some between-subject data that could be used for such tests. However, no statistics regarding between-subject results are explicitly reported in those papers. This may in part be due to the higher number of subjects needed in between-subject designs to reach a level of statistical power that is comparable to within-subject designs. Indeed, by removing subject variance from the error term, within-subject tests are statistically much more powerful than between-subject designs. This means that not only a between-subject design will need n x k the subjects of a within-subject design with k subjects ( where n is the number of treatments), but the within-subject design using k subjects will still be much more powerful [START_REF] Greenwald | Within-Subject Designs: To Use or Not to Use[END_REF]. This may in part explain the recurrence to within-subject designs, inasmuch as obtaining between-subject data on the effect of high stakes can quickly become very expensive. Even though we do not agree with certain methodological arguments that consider between-subject data the gold standard to which all other results need to be compared (see e.g. [START_REF] Poulton | Unwanted Range Effects from Using Within-Subject Experimental Designs[END_REF], we think that such data may in some cases reveal interesting new insights.

We want to add to the existing literature by presenting some between-subject results on the effect of financial stakes on risk attitudes for small probability prospects. Even though our attention is restricted to small probabilities for budgetary reasons, the changes in stakes are substantial, ranging from prizes of €4 ($6) to prizes of €100 ($150). Also, we collected data using different methodologies, ranging from certainty equivalents to buying prices (amounts to invest). While finding the typical pattern of risk seeking for small probability prospects under low monetary stakes, we show that such risk seeking is substantially reduced under high stakes.

THE EXPERIMENT

Method

Subjects. The experiment has been conducted at GATE (Groupe d'Analyse et de Théorie Economique), at the University of Lyon, France. 78 undergraduate students from the local engineering and business schools were recruited using the ORSEE software [START_REF] Greiner | An Online Recruitment System for Economic Experiments[END_REF]. 64% of subjects were female, the average age was 22. Four sessions were run, with 20 subjects in two of them and 19 subjects in the other two. On average the subjects earned €22.58, including a show-up fee of €5, for an experiment lasting less than 30 minutes.

Tasks. The experiment was computerized, using the REGATE software [START_REF] Zeiliger | A Presentation of Regate, Internet Based Software for Experimental Economics[END_REF].

Subjects absolved several tasks in the course of an experiment on lottery preferences and probability representations that is described in [START_REF] Lefebvre | The Ratio Bias Phenomenon: Fact or Artifact? Working Paper[END_REF]. Since different probability representations were found to have no significant effect, we will not further discuss how probabilities were represented. Two different tasks relevant for the present paper were used. First, a context-free neutral task was used to elicit certainty equivalents for a binary lottery giving a 10% probability of winning a prize and a complementary probability of winning nothing. Certainty equivalents were elicited through a list of 26 choices (see Appendix).

Next, subjects' willingness to invest into a risky and potentially lucrative project was explored. The investment task was presented as the decision to invest or not in a clean energy project. Subjects were given an initial endowment, and the amount subjects were willing to invest was then elicited through a list of 12 choices. In the investment task, subjects were randomly assigned one of three probabilities of investment success ranging between 5.9% and 7% . Since this was done for both the low and high stakes conditions, and since ratios of elicited values to expected value (EV) are used for the analysis, this small variation in probabilities does not affect the results presented.

Incentives. A show-up fee of €5 was provided to all subjects. Subjects were assigned to either a Low-Stakes or a High-Stakes condition. In the Low-Stakes condition, the prizes were €10 ($15) for the neutral task, and €4 ($6) for the investment task, with the latter to be financed out of an initial endowment of €0.60 (90¢). In the High-Stakes condition, all amounts were increased by a factor of 10, implying prizes of €100 ($150) and €40 ($60) for the neutral and investment tasks, respectively.

Encoding. The certainty equivalent (CE) and the willingness-to-pay for investment (WTP) were calculated as the mean between the two amounts for which subjects switched from the prospect to the certain amount (in the case of CE), or from the certain amount to the prospect (in the case of WTP).

Hypotheses. Given that small probabilities are used throughout, we expect that subjects will on average be risk seeking as predicted by prospect theory [START_REF] Abdellaoui | Parameter-Free Elicitation of Utility and Probability Weighting Functions[END_REF][START_REF] Bleichrodt | A Parameter-Free Elicitation of the Probability Weighting Function in Medical Decision Analysis[END_REF][START_REF] Van De Kuilen | A Midpoint Technique for Easily Measuring Prospect Theory's Probability Weighting[END_REF][START_REF] Wu | Curvature of the Probability Weighting Function[END_REF]. Indeed, subjects are generally found to overweight small probabilities of winning, so that we may expect CEs and WTPs to lie above the expected value. Finally and most importantly, we expect that risk seeking will be reduced significantly in the High-Stakes condition as compared to the Low-Stakes condition.

Results

Neutral task. 6 subjects were dropped from the sample because they switched multiple times between the sure amount and the prospect. Those subjects are dropped since it is not clear what their CE is. Under low stakes, we find the typical behavior of risk seeking for small probabilities predicted by prospect theory, with a mean ratio of the CE to EV of 1.66. We thus strongly reject the hypothesis that subjects are expected value maximizers for low stakes (t(39) = 8.53, p < 0.001; all p-values are two-sided) in favor of the hypothesis that subjects are risk seeking. At the individual level, only two subjects out of 40 can be classified as risk averse, three as risk neutral, and the remaining 35 as risk seeking. Data on the frequency distribution of risk attitudes are displayed in Table 1. This picture changes when we look at the High-Stakes condition. On average, subject are now roughly risk neutral with an average ratio of CE to EV of 0.96. Indeed, the hypothesis of risk neutrality can now not be rejected (t(31) = -0.44, p = 0.67). At the individual level, 14

subjects can now be classified as risk averse, six as risk neutral, and 12 as risk seeking. We thus confirm that in the High-Stakes condition subjects are on average significantly less risk seeking than in the Low-Stakes condition (z = 4.74, p < 0.001, Mann-Whitney test). Though women are on average slightly less risk seeking than men, this difference is neither significant for the overall data nor within either of the two incentive conditions. Again, we find the typical behavior of risk seeking for low stakes, with a mean ratio of WTP to EV equal to 1.88. We again easily reject the null hypothesis of risk neutrality (t(37) = 6.05, p < 0.001). At the individual level, five subjects can be classified as risk averse, two as risk neutral, and 33 as risk seeking (see Table 1). Under high stakes we now also find risk seeking, with the mean WTP on EV ratio equal to 1.52. This time we reject the null hypothesis of risk neutrality also for high stakes (t(33) = 3.31, p = 0.002). At the individual level, we can classify nine subjects as risk averse, two as risk neutral, and 23 as risk seeking. As hypothesized, subjects in the High-Stakes condition are on average less risk seeking than subjects in the Low-Stakes condition (z = 1.93, p = 0.05; Mann-Whitney test) (see also Figure 2). Just as for the neutral task, no significant gender differences are present.

Figure 2: Proportion of subjects by WTP/EV ratio. Approximations of WTP/EV were used to condensate findings; this was necessary since different probability levels produced small difference in EV that resulted in a multiplicity of slightly different WTP/EV ratios.

DISCUSSION

Our between-subject data confirm previous findings in the literature according to which individuals become more risk averse-or in our case less risk seeking-when high stakes are involved. This effect is indeed very strong for a neutral task in which CEs are elicited. When smaller stakes are used in an investment task, the effect is found also for willingness to invest into a potentially profitable project. However, the effect is significantly less strong in the latter case (Z = 1.49, p = 0.06, Fisher's z test).

Subjects are generally found to be risk seeking for small probability prospects as would be predicted by the overweighting of small probabilities generally found in the literature [START_REF] Abdellaoui | Parameter-Free Elicitation of Utility and Probability Weighting Functions[END_REF][START_REF] Bleichrodt | A Parameter-Free Elicitation of the Probability Weighting Function in Medical Decision Analysis[END_REF][START_REF] Van De Kuilen | A Midpoint Technique for Easily Measuring Prospect Theory's Probability Weighting[END_REF][START_REF] Wu | Curvature of the Probability Weighting Function[END_REF]. While subjects become risk neutral under high incentives for the neutral task, they remain risk seeking for the investment task. This difference between the neutral and the investment tasks may potentially be driven by either of two factors: the introduction of context into the decision (people may be more prone to risk-seeking in a context of investment in an environmental good), or the elicitation of willingness-to-pay instead of certain equivalents. Disentangling these two explanations is beyond the scope of this paper but may be an interesting topic for future research.

Assuming prospect theory as a descriptive model of choice, the difference between low and high stakes can be explained with attitudes towards outcomes, since probability weighting is a purely probabilistic matter and should not change between low and high stakes (see however Bradley, 2003, on difficulties with this separation). This is however somewhat troubling, since it is generally assumed that utility should be linear for such small amounts [START_REF] Abdellaoui | A Tractable Method to Measure Utility and Loss Aversion under Prospect Theory[END_REF][START_REF] Booij | A Parameter-Free Analysis of the Utility of Money for the General Population under Prospect Theory[END_REF]. One potential explanation is that the low income of our subject population of students implies utility curvature already for relatively small gains. A potential alternative explanation would be that utility is not consistently concave over the outcome space, but rather contains a convex section for very small amounts [START_REF] Bosch-Domènech | Does Risk Aversion or Attraction Depend on Income? An Experiment[END_REF][START_REF] Markowitz | The Utility of Wealth[END_REF].

According to this peanut effect [START_REF] Prelec | Decision Making over Time and Under Uncertainty: A Common Approach[END_REF][START_REF] Weber | Playing for Peanuts: Why is risk-seeking more common for low-stakes games? Organizational Behavior and Human Decision Processes[END_REF], the extreme risk seeking found for small amounts would be the result of both overweighing of small probabilities and an increasing marginal utility of money for small amounts of money. Indeed, the latter would predict that small amounts of money that are offered for sure in our two tasks are attributed very low utility, which implies extreme risk seeking in the Low-Stakes condition. This approach would explain the extreme difference in risk seeking behavior in our neutral task not through a rapidly decreasing marginal utility of money, but rather through a combination of initially increasing and subsequently decreasing marginal utility (in addition to the overweighting of small probabilities). Notice how this conceptual framework may also explain why we find a weaker effect of monetary stakes in 1 0 the investment task. Indeed, sure amounts used in the latter are even smaller and may thus be undervalued in both conditions, resulting in persisting risk seeking in the High-Stakes condition.

CONCLUSION

The importance of monetary stakes for risk attitude is an important issue in the decision making literature. Even though several studies investigating the effect of high monetary stakes on risk attitude exist, data reported are generally obtained by means of within-subject designs, and doubts have been aired about the soundness of that approach. In this paper we test instead the effect of high monetary stakes for small probability prospects in a between-subject design.

While finding the typical pattern of risk seeking for small probability prospects under low stakes, such risk seeking is found to be substantially reduced under high stakes. This finding may create problems for prospect theory as it stands, but could be explained either by subject pool effects, or by a peanut effect, according to which the utility curve may have a convex section for low gain amounts.

APPENDIX

Neutral Task (High stakes):

Figure 1 :

 1 Figure 1: Ratio of CE to EV for Low and High Stakes

Table 1 :

 1 Classification of subjects in terms of risk attitude

	Tasks	Condition	Distribution of risk attitudes
			Risk	Risk	Risk	Total
			averse	neutral	seeking
	Neutral task	Low stakes	2 (5.00)	3 (7.50)	35 (87.50)	40 (100)
		High stakes	14 (43.75)	6 (18.75)	12 (37.50)	32 (100)
	Investment task	Low stakes	5 (12.50)	2 (5.00)	33 (82.50)	40 (100)
		High stakes	9 (26.47)	2 (5.88)	23 (67.65)	34 (100)
	Note: Relative frequencies in parentheses.

The authors thank S. Ferriol for programming the experiment. Financial support from the EMIR program of the French National Agency for Research (ANR) is gratefully acknowledged (BLAN07-3_185547).