Andrei Zaharescu

Radu Horaud
email: radu.horaud@inrialpes.fr

Robust Factorization Methods Using A Gaussian/Uniform Mixture Model

Keywords: robust factorization, 3-D reconstruction, multiple camera calibration, data clustering, expectation-maximization, EM, M-estimators, outlier rejection

In this paper we address the problem of building a class of robust factorization algorithms that solve for the shape and motion parameters with both affine (weak perspective) and perspective camera models. We introduce a Gaussian/uniform mixture model and its associated EM algorithm. This allows us to address robust parameter estimation within a data clustering approach. We propose a robust technique that works with any affine factorization method and makes it robust to outliers. In addition, we show how such a framework can be further embedded into an iterative perspective factorization scheme. We carry out a large number of experiments to validate our algorithms and to compare them with existing ones. We also compare our approach with factorization methods that use M-estimators.

Introduction

The problem of 3-D reconstruction from multiple images is central in computer vision [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF][START_REF] Luong | The Geometry of Multiple Images[END_REF].

Bundle adjustment provides both a general method and practical algorithms for solving this reconstruction problem using maximum likelihood [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]. Nevertheless, bundle adjustment is nonlinear in nature and sophisticated optimization techniques are necessary, which in turn require proper initialization. Moreover, the combination of bundle adjustment with robust statistical methods to reject outliers is not clear both from the points of view of convergence properties and of efficiency. Factorization was introduced by Tomasi & Kanade [START_REF] Tomasi | Shape and motion from image streams under orthography: a factorization method[END_REF] as an elegant solution to affine multiple-view reconstruction; their initial solution based on SVD and on a weakperspective camera model has subsequently been improved and elaborated by Morris & Kanade [START_REF] Morris | A unified factorization algorithm for points, line segments and planes with uncertainty models[END_REF], Anandan & Irani [START_REF] Anandan | Factorization with uncertainty[END_REF], Hartley & Schaffalitzky [START_REF] Hartley | Powerfactorization: 3d reconstruction with missing or uncertain data[END_REF] as well as by many others. These methods treat the non-degenarate cases. Kanatani [START_REF] Kanatani | Geometric information criterion for model selection[END_REF], [START_REF] Kanatani | Uncertainty modeling and model selection for geometric inference[END_REF] investigated how to apply model selection techniques to deal with degenarate cases, namely when the 3-D points lie on a plane and/or the camera centers lie on a circle.

The problem can be formulated as the one of minimizing the following Frobenius norm:

θ * = arg min θ S -Ŝ(θ) 2 F ( 1
)
where matrix S = [s ij ] denotes the measurement matrix containing matched 2-D image observations, Ŝ(θ) = MP denotes the prediction matrix that can be factorized into the affine motion matrix M and the affine shape matrix P. Hence, we denote by θ the affine motion and shape parameters collectively. In the error-free case, direct factorization of the observation matrix using SVD provides an optimal solution. More recently the problem of robust affine factorization has received a lot of attention and powerful algorithms that can deal with noisy, missing, and/or erroneous data were suggested.

Anandan & Irani [START_REF] Anandan | Factorization with uncertainty[END_REF] extended the classical SVD approach to deal with the case of directional uncertainty. They used the Mahalanobis norm instead of the Frobenius norm and they reformulated the factorization problem such that the Mahalanobis norm can be transformed into a Frobenius norm. This algorithm handles image observations with covariance up to a few pixels but it cannot cope with missing data, mismatched points, and/or outliers. More generally, a central idea is to introduce a weight matrix W of the same size as the measurement matrix S. The minimization criterion then becomes:

θ * = arg min θ W ⊗ (S -Ŝ(θ)) 2 F (2)
where ⊗ denotes the Hadamard product (A = B ⊗ C ⇐⇒ a ij = b ij c ij ) and W = [w ij ] is matrix whose entries are weights that reflect the confidence associated with each image observation. The most common way of minimizing eq. ( 2) is to use alternation methods: these methods are based on the fact that, if either one of the matrices M or P is known, then there is a closed-form solution for the other matrix that minimizes eq. ( 2). Morris & Kanade [START_REF] Morris | A unified factorization algorithm for points, line segments and planes with uncertainty models[END_REF] were the first to propose such an alternation method. The PowerFactorization method introduced by Hartley & Schaffalitzky [START_REF] Hartley | Powerfactorization: 3d reconstruction with missing or uncertain data[END_REF], as well as other methods by Vidal & Hartley [START_REF] Vidal | Motion segmentation with missing data using powerfactorization and gpca[END_REF], and Brant [START_REF] Brant | Closed-form solutions for affine reconstruction under missing data[END_REF] fall into this category.

PowerFactorization is based on the PowerMethod for sparse matrix decomposition [START_REF] Golub | Matrix Computations[END_REF]. Notice that these techniques are very similar in spirit with PCA methods with missing data, Wiberg [START_REF] Wiberg | Computation of principal components when data are missing[END_REF],

Ikeuchi, Shum, & Reddy [START_REF] Shum | Principal component analysis with missing data and its application to polyhedral object modeling[END_REF], Roweis [START_REF] Roweis | EM algorithm for PCA and SPCA[END_REF], and Bishop [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. Another way to alternate between the estimation of motion and shape is to use factor analysis, Gruber and Weiss [START_REF] Gruber | Factorization with uncertainty and missing data: exploring temporal coherence[END_REF], [START_REF] Gruber | Multibody factorization with uncertainty and missing data using the em algorithm[END_REF].

Alternatively, robustness may be achieved through adaptive weighting, i.e., by iteratively updating the weight matrix W which amounts to modifying the data S, such as is done by Aanaes et al. [START_REF] Aanaes | Robust factorization[END_REF]. Their method uses eq. ( 1) in conjunction with a robust loss function (see Stewart [START_REF] Stewart | Robust parameter estimation in computer vision[END_REF] and Meer [START_REF] Meer | Robust techniques for computer vision[END_REF] for details) to iteratively approximate eq. ( 2), getting a temporary optimum. The approximation is performed by modifying the original data S such that the solution to eq. [START_REF] Aanaes | Robust factorization[END_REF] with modified data S is the same as the solution to eq. ( 2) with the original data:

θ * = arg min θ W ⊗ (S -Ŝ(θ)) 2 F = arg min θ S -Ŝ(θ)) 2 F (3)
In [START_REF] Aanaes | Robust factorization[END_REF] the weights are updated via IRLS [START_REF] Stewart | Robust parameter estimation in computer vision[END_REF]. This may well be viewed as both an iterative and an alternation method because the motion matrix M is estimated using SVD, then the shape matrix P is estimated knowing M, while the image residuals are calculated and the data (the weights) are modified, etc. A similar algorithm that performs outlier correction was proposed by Huynh, Hartley, & Heyden [START_REF] Huynh | Outlier correction in image sequences for the affine camera[END_REF]. Indeed, if the observations are noisy, the influence of outliers can be decreased by iteratively replacing bad observations with "pseudo" observations. The convergence of such methods, as [START_REF] Aanaes | Robust factorization[END_REF] or [START_REF] Huynh | Outlier correction in image sequences for the affine camera[END_REF] is not proved but is tested through experiments with both simulated and real data.

An alternative to M-estimators are random sampling techniques developed independently in computer vision [START_REF] Fischler | Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography[END_REF] and statistics [START_REF] Rousseeuw | Least median of squares regression[END_REF] (see Meer [START_REF] Meer | Robust techniques for computer vision[END_REF] for a recent overview of these methods).

For example, Huynh & Heyden [START_REF] Huynh | Robust factorization for the affine camera: Analysis and comparison[END_REF] and Tardif et al. [START_REF] Ph | Algorithms for batch matrix factorization with application to structure-from-motion[END_REF] use RANSAC, Trajkovic and Hedley use LMedS [START_REF] Trajkovic | Robust recursive structure and motion recovery under affine projection[END_REF], and Hajder and Chetverikov [START_REF] Hajder | Robust structure from motion under weak perspective[END_REF] use LTS (Least Trimmed Squares) [START_REF] Rousseeuw | Positive-breakdown robust methods in computer vision[END_REF]. The major drawback of these methods is that they must consider a large number of subsets sampled from the observation matrix S.

Generally speaking, robust regression techniques, such as the ones that we briefly discussed, work well in conjunction with affine factorization algorithms. Factorization was initially designed as a "closed-form solution" to multiple-view reconstruction, but robust affine factorization methods are iterative in nature, as explained above. This has several implications and some drawbacks. In the presence of a large number of outliers, proper initialization is required. The use of an influence function (such as the truncated quadratic) that tends to zero too quickly cause outliers to be ignored and hence, this raises the question of a proper choice of an influence function. The objective function is non-convex implying that IRLS will be trapped in local minima.

The generalization of affine factorization to deal with perspective implies the estimation of depth values associated with each reconstructed point. This is generally performed iteratively [START_REF] Sturm | A factorization based algorithm for multi-image projective structure and motion[END_REF], [START_REF] Christy | Euclidean shape and motion from multiple perspective views by affine iterations[END_REF], [START_REF] Mahamud | Iterative projective reconstruction from multiple views[END_REF], [START_REF] Mahamud | Provably-convergent iterative methods for projective structure from motion[END_REF], [START_REF] Miyagawa | Motion and shape recovery based on iterative stabilization for modest deviation from planar motion[END_REF], [START_REF] Oliensis | Iterative extensions of the sturm/triggs algorithm: Convergence and nonconvergence[END_REF]. It is not yet clear at all how to combine iterative robust methods with iterative projective/perspective factorization methods.

In this paper we cast the problem of robust factorization into the framework of data clustering [START_REF] Fraley | Model-based clustering, discriminant analysis, and density estimation[END_REF]. Namely, we consider the problem of classifying the observed 2-D matched points into two categories: inliers and outliers. For that purpose we model the likelihood of the observations with a Gaussian/uniform mixture model. This leads to a maximum likelihood formulation with missing variables that can be solved with the EM algorithm [START_REF] Dempster | Maximum likelihood estimation from incomplete data via the EM algorithm (with discussion)[END_REF], [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF], [START_REF] Fraley | Model-based clustering, discriminant analysis, and density estimation[END_REF]. Notice that this approach is different than the method proposed by Miller & Browning [START_REF] Miller | A mixture model and em-based algorithm for class discovery, robust classification and outlier rejection in mixed labeled/unlabeled data sets[END_REF] requiring both labeled and ulabeled data sets.

We devise an EM algorithm within the framework of 3-D reconstruction and within the specific mixture model just outlined. Each EM step guarantees that the likelihood is increased.

Hence EM may indeed lead to a local maximum of the likelihood. We show that in this particular case (normally distributed inliers and uniformly distributed outliers) the posterior propabilities have a very simple interpretation in terms of robust regression. We describe an affine factorization algorithm that uses EM; This algorithm is robust and it shares the convergence properties just outlined. We also describe an extension of this algorithm to deal with the perspective camera model.

We performed several experiments in two different scenarios: multiple-camera calibration and 3-D reconstruction using turn-table data. Our method was compared to other methods on an equal footing: it performs as well as bundle adjustment to estimate external camera parameters. It performs better than IRLS (used in conjunction with the truncated quadratic) to eliminate outliers in some difficult cases.

The remainder of this paper is organized as follows. Section 2 describes the probabilistic modelling of inliers and outliers using a mixture between a Gaussian and an uniform distribution.

Section 3 explains how this probabilistic model can be used to derive an affine factorization algorithm and section 4 extends this algorithm to iterative perspective factorization. Sections 5 and 6 describe experiments performed with multiple-camera calibration and with multi-view reconstruction data sets. Section 7 compares our approach to M-estimators and section 8 draws some conclusions and gives some directions for future work.

Probabilistic modelling of inlier/outlier detection

The 2-D image points s ij (1 ≤ i ≤ k, 1 ≤ j ≤ n) are the observed values of an equal number of random variables s ij . We introduce another set of random variables, z ij which assign a category to each observation. Namely there are two possible categories, an inlier category and an outlier category. More specifically z ij = inlier means that the observation s ij is an inlier while z ij = outlier means that the observation s ij is an outlier.

We define the prior probabilities as follows. Let A i be the area associated with image i and we assume that all the images have the same area, A i = A. If an observation is an inlier, then it is expected to lie within a small circular image patch a of radius σ 0 , a = πσ 2 0 . The prior probability of an inlier is the proportion of the image restricted to such a small circular patch:

P (z ij = inlier) = a A (4)
Similarly, if the observation is an outlier, its probability should describe the fact that it lies outside this small patch:

P (z ij = outlier) = A -a A (5)
Moreover, an observation s ij , given that it is an inlier, should lie in the neighborhood of an estimation ŝij . Therefore, we will model the probability of an observation s ij given that it is assigned to the inlier category with a Gaussian distribution centered on ŝij and with a 2×2 covariance matrix C. We obtain:

P θ (s ij |z ij = inlier) = 1 2π(det C) 1/2 exp - 1 2 d 2 C (s ij , ŝij (θ)) (6)
where we denote by d C the Mahalanobis distance:

d 2 C (s ij , ŝij (θ)) = (s ij -ŝij (θ)) ⊤ C -1 (s ij -ŝij (θ)) (7)
Whenever the observation is an outlier, it may lie anywhere in the image. Therefore, we will model the probability of an observation s ij given that it is assigned to the outlier category with a uniform distribution over the image area:

P (s ij |z ij = outlier) = 1 A (8)
Since each variable z ij can take only two values, marginalization is straightforward and we obtain:

P θ (s ij ) = P θ (s ij |z ij = inlier)P (z ij = inlier) + P (s ij |z ij = outlier)P (z ij = outlier) = a 2π(det C) 1/2 A exp - 1 2 d 2 C (s ij , ŝij (θ)) + A -a A 2 (9)
We already defined the small area a as a disk of radius σ 0 , a = πσ 2 0 and we assume that a ≪ A. Using Bayes' formula 1 , we obtain the posterior conditional probability of an observation to be an inlier. We denote this posterior probability by α in ij :

α in ij = P θ (z ij = inlier|s ij ) = 1 1 + 2 σ 2 0 (det C) 1/2 exp 1 2 d 2 C (s ij , ŝij (θ)) (10)
The covariance matrix can be written as C = U ΛU ⊤ where U is a rotation and Λ is a diagonal form with entries λ 1 and λ 2 . Hence det(C) = λ 1 λ 2 . In order to plot and illustrate the shape of α in ij as a function of C we consider the case of an isotropic covariance, i.e., λ 1 = λ 2 = σ 2 and one may notice that the rotation becomes irrelevant in this case. We have: C = σ 2 I 2 . Eq. [START_REF] Fraley | Model-based clustering, discriminant analysis, and density estimation[END_REF] writes in this case:

α in ij = P θ (z ij = inlier|s ij ) = 1 1 + 2σ 2 σ 2 0 exp s ij -ŝij (θ) 2 2σ 2 (11)
This posterior probability is shown on Figure 1, i.e., the function f σ (x) = 1/(1+σ 2 exp(x 2 /2σ 2 )).

Here σ takes discrete values in the interval [0.05, 5] and σ 2 0 = 2, i.e., inliers lie within a circle of radius 2 pixels centered on a prediction. It is worthwhile to notice that, at the limit σ → 0, we obtain a Dirac function:

f 0 (x) =    1 if x = 0 0 if x = 0 (12)
Figure 1: Various plots of the conditional posterior probability of an observation to be an inlier, i.e., f σ (x) = 1/(1 + σ 2 exp(x 2 /2σ 2 )). This function corresponds to eq. ( 11) with σ 2 0 = 2. As the variance decreases, i.e., σ = 5, 4, 3, 2, 1, 0.5, 0.25, 0.1, 0.05, the function becomes more and more discriminant. It is worthwhile to notice that lim σ→0 f σ (x) is a Dirac.

The posterior conditional probability of an observation to be an outlier is given by:

α out ij = P θ (z ij = outlier|s ij ) = 1 -α in ij (13)

Maximum likelihood with inliers

The maximum likelihood estimator (ML) maximizes the log-likelihood of the joint probability of the set of measurements, P θ (S). Under the assumption that the observations are independent and identically distributed we have:

P θ (S) = i,j P θ (s ij ) (14)
Since we assume that all the observations are inliers, eq. ( 9) reduces to:

P θ (s ij ) = P θ (s ij |z ij = inlier) (15)
The log-likelihood of the joint probability becomes:

log P θ (S) = - 1 2 i,j d 2 C (s ij , ŝij (θ)) + log(det C) + const (16)
which can be written as the following criterion:

Q M L = 1 2 i,j (s ij -ŝij (θ)) ⊤ C -1 (s ij -ŝij (θ)) + log(det C) (17)
The shape and motion parameters can be estimated by minimizing the above criterion with respect to θ:

θ * = arg min θ 1 2 i,j (s ij -ŝij (θ)) ⊤ C -1 (s ij -ŝij (θ)) (18)
Once an optimal solution is found, i.e., θ * , it is possible to minimize eq. ( 17) with respect to the covariance matrix which yields (see appendix A):

C * = 1 m i,j (s ij -ŝij (θ * ))(s ij -ŝij (θ * )) ⊤ (19)
where m = k × n is the total number of observations for k images and n 3-D points.

Alternatively, if one uses an isotropic covariance, i.e., C = σ 2 I, By minimization of Q M L with respect to θ we obtain:

θ * = arg min θ 1 2 i,j s ij -ŝij (θ)) 2 (20)
The optimal variance is given by (see appendix B):

σ 2 * = 1 2m i,j s ij -ŝij (θ * )) 2 (21)

Maximum likelihood with inliers and outliers

In the presence of outliers, the previous method cannot be applied. Instead, one has to use the joint probability of the observations and of their assignments. Again, by assuming that the observations are independent, we have:

P θ (S, Z) = i,j P θ (s ij , z ij ) = i,j P θ (s ij |z ij )P (z ij ) = i,j P θ (s ij |z ij = inlier)P (z ij = inlier) δ in (z ij ) P (s ij |z ij = outlier)P (z ij = outlier) δout(z ij ) (22)
The random variables δ in (z ij ) and δ out (z ij ) are defined by:

δ in (z ij ) =    1 if z ij = inlier 0 otherwise δ out (z ij ) =    1 if z ij = outlier 0 otherwise
By taking the logarithm of the above expression and grouping constant terms, we obtain:

log P θ (S, Z) = i,j δ in (z ij ) log(P θ (s ij |z ij = inlier)) + δ out (z ij ) log(P θ (s ij |z ij = outlier)) + const (23)
This cannot be solved as previously becasue of the presence of the missing assignment variables z ij . Therefore, they will be treated within an expectation-maximization framework. For this purpose we evaluate the conditional expectation of the log-likelihood over the random variables z ij , given the observations S:

E Z log(P θ (S, Z))|S = i,j log(P θ (s ij |z ij = inlier))E Z [δ in (z ij )|S] + log(P (s ij |z ij = outlier))E Z [δ out (z ij )|S] (24)
In this formula we omitted the constant terms, i.e., the terms that do not depend on the parameters θ and C. The subscript Z indicates that the expectation is taken over the random variable z.

From the definition of δ in (z ij ) we have:

E[δ in (z ij )] = δ in (z ij = inlier)P (z ij = inlier) + δ in (z ij = outlier)P (z ij = outlier) = P (z ij = inlier)
Hence:

E Z [δ in (z ij )|S] = P (z ij = inlier|S) = P (z ij = inlier|s ij ) = α in ij
and:

E Z [δ out (z ij )|S] = 1 -α in ij
Therefore, after removing constant terms, the conditional expectation becomes:

E Z log(P θ (S, Z))|S = - 1 2 i,j α in ij d 2 C (s ij , ŝij (θ)) + log(det C) (25)
This leads to the following criterion:

Q EM = 1 2 i,j α in ij (s ij -ŝij (θ)) ⊤ C -1 (s ij -ŝij (θ) + log(det C) (26)

Robust affine factorization with the EM algorithm

In this section we provide the details of how the robust affine factorization problem can be solved iteratively by maximum likelihood via the expectation-maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood estimation from incomplete data via the EM algorithm (with discussion)[END_REF], and how the observations can be classified into either inliers or outliers by maximum a-posteriori (MAP).

By inspection of equations ( 17) and ( 26) one may observe that the latter is a weighted version of the former and hence our formulation has strong simililarities with M-estimators and their practical solution, namely iteratively reweighted least-squres (IRLS) [START_REF] Stewart | Robust parameter estimation in computer vision[END_REF]. Nevertheless, the weights ω ij = α in ij were obtained using a Bayesian approach: they correspond to the posterior conditional probabilities of the observations (i.e., given that they are inliers), and such that the equality α in ij +α out ij = 1 holds for each observation. The structure and the shape of these posteriors are depicted by equations ( 10) and ( 11) and shown on Figure 1. These probabilities are functions of the residual but they are parameterized as well by the 2×2 covariance matrix C associated with the normal probability distribution of the observations: One advantage of our formulation over IRLS is that this covariance is explicitly taken into consideration and estimated within the EM algorithm.

It is worthwhile to remark that the minimization of eq. ( 18) over the affine shape and motion parameters, i.e., θ, can be solved using an affine camera model and a factorization method such that the ones proposed in the literature [START_REF] Aanaes | Robust factorization[END_REF][START_REF] Hartley | Powerfactorization: 3d reconstruction with missing or uncertain data[END_REF]. In practice we use the PowerFactorization method proposed in [START_REF] Hartley | Powerfactorization: 3d reconstruction with missing or uncertain data[END_REF]. The minimization of eq. ( 26) can be solved in the same way, provided that estimates for the posterior probabilities α in ij are available. This can be done by iterations of the EM algorithm:

• The E-step computes the conditional expectation over the assignment variables associated with each observation, i.e., eq. ( 25). This requires a current estimate of both θ and C from which the α in ij 's are updated.

• The M-step maximizes the conditional expectation or, equivalently, minimizes eq. ( 26)

with fixed posterior probabilities. This is analogous, but not identical, with finding the means µ ij and a common covariance

C of m = k × n Gaussian distributions, with µ ij = ŝij (θ).
Nevertheless, the means µ = {µ 11 , . . . , µ kn } are parameterized by the global variables θ. For this reason, the minimization problem needs a specific treatment (unlike the classical mixture of Gaussians approach where the means are independent).

Therefore min µ Q EM in the standard EM method must be replaced by min θ Q EM and it does depend on C in this case:

θ * = arg min θ 1 2 i,j α in ij (s ij -ŝij (θ)) ⊤ C -1 (s ij -ŝij (θ)) (27)
Moreover, the covariance that minimizes eq. ( 26) can be easily derived from [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]:

C * = 1 i,j α in ij i,j α in ij (s ij -ŝij (θ * ))(s ij -ŝij (θ * )) ⊤ (28)
In many practical situations it is worthwhile to consider the case of an isotropic covariance, in which case the equations above reduce to:

θ * = arg min θ 1 2 i,j α in ij s ij -ŝij (θ)) 2 (29)
and

σ 2 * = 1 2 i,j α in ij i,j α in ij s ij -ŝij (θ * )) 2 (30)
This may well be viewed as a special case of model-based clustering [START_REF] Fraley | Model-based clustering, discriminant analysis, and density estimation[END_REF]. It was proved [START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF] that EM guarantees convergence, i.e., that Q q+1 EM < Q q EM , where the overscript q denotes the q th iteration, and that this implies the maximization of the joint probability of the observations: P θ (S) q+1 > P θ (S) q . To conclude, the algorithm can be paraphrased as follows:

Affine factorization with EM:

Initialization: Use the PowerFactorization method to minimize eq. ( 18). This provides initial estimates for θ (the affine shape and motion parameters). Estimate C (the covariance matrix) using eq. [START_REF] Huynh | Robust factorization for the affine camera: Analysis and comparison[END_REF].

Iterate until convergence

Expectation: Update the values of α in ij according to eq. ( 10) or eq. ( 11).

Maximization: Minimize Q EM over θ (affine factorization) using either eq. ( 27)

or eq. ( 29). Compute the covariance C with eq. ( 28) or the variance σ 2 with eq. [START_REF] Morris | A unified factorization algorithm for points, line segments and planes with uncertainty models[END_REF].

Maximum a posteriori: Once the EM iterations terminate, choose in between inlier and outlier for each observation, i.e., max{α in ij ; α out ij }.

The algorithm needs initial estimates for the shape and motion parameters from which an initial covariance matrix can be estimated. This guarantees that, at the start of EM, all the residuals have equal importance. Nevertheless, "bad" observations will have a large associated residual and, consequently, the covariance is proportionally large. As the algorithm proceeds, the covariance adjusts to the current solution while the posterior probabilities α in ij become more and more discriminant as depticted on Figure 1. Eventually, observations associated with small residuals will be classified as inliers, and observations with large residuals will be classified as outliers.

The overall goal of 3-D reconstruction consists of the estimation of the shape and motion parameters: As just explained, we embed affine reconstruction in the M-step. Therefore, with our algorithm, robustness stays outside the factorization method at hand -is it iterative or notand hence one can plug into EM any factorization procedure.

Robust perspective factorization

In this section we address the problem of 3-D reconstruction using instrinsically calibrated cameras. Moreover, we consider both the weak-perspective and the perspective camera models, and we explain how the affine solution provided by factorization can be upgraded to Euclidean reconstruction. We describe an algorithm that combines the EM affine factorization algorithm described above with an iterative perspective factorization algorithm [START_REF] Christy | Euclidean shape and motion from multiple perspective views by affine iterations[END_REF][START_REF] Zaharescu | Multiple camera calibration using robust perspective factorization[END_REF]. This results in a robust method for solving the 3-D Euclidean reconstruction problem as well as the multiple-camera calibration problem. An image point s = (x, y) is the projection of a 3-D point X:

x ij = r x i • Xj + t x i r z i • Xj + t z i = a x i • Xj + b x i ε ij + 1 ( 31
)
y ij = r y i • Xj + t y i r z i • Xj + t z i = a y i • Xj + b y i ε ij + 1 (32)
We introduced the following notations: The rotation matrix R ⊤ i = [r x i r y i r z i ] and the translation vector t ⊤ i = (t x i t y i t z i ) correspond to the motion parameters and they are also denoted the external camera parameters. Dividing the above equations with the depth t z i we obtain a similar set of scaled equations. We have:

a x i = r x i /t z i , a y i = r y i /t z i , b x i = t x i /t z i and b y i = t y i /t z i .
We denote by ε ij the perspective distorsion parameters, namely the following ratios:

ε ij = r z i • Xj t z i (33)
Finally, the perspective equations, i.e., eqs. ( 31) and ( 32) can be written as:

s ij (1 + ε ij ) = A i X j (34)
where X = ( X, 1) and A i denotes the following 2×4 matrix:

A i =   a x i b x i a y i b y i  
From now on we can replace the parameter vector θ with the affine shape and motion parameters, namely the point set X = {X 1 , . . . X j , . . . X k } and the matrix set

A = {A 1 , . . . A j , . . . A n }.
Using these notations, eq. ( 27) can now be written as:

min A,X 1 2 i,j α in ij (s ij (1 + ε ij ) -A i X j ) ⊤ C -1 (s ij (1 + ε ij ) -A i X j ) (35)
which can be solved via the EM affine factorization algorithm with ε ij = 0, ∀(i, j). A weakperspective camera model can then be used for upgrading to Euclidean reconstruction.

The introduction of the perspective camera model adds non null perspective-distorsion parameters ε ij , i.e., eq. ( 33). One fundamental observation is the following: If estimates for the both with uncalibrated and calibrated cameras [START_REF] Sturm | A factorization based algorithm for multi-image projective structure and motion[END_REF], [START_REF] Christy | Euclidean shape and motion from multiple perspective views by affine iterations[END_REF], [START_REF] Mahamud | Iterative projective reconstruction from multiple views[END_REF], [START_REF] Mahamud | Provably-convergent iterative methods for projective structure from motion[END_REF], [START_REF] Miyagawa | Motion and shape recovery based on iterative stabilization for modest deviation from planar motion[END_REF], [START_REF] Oliensis | Iterative extensions of the sturm/triggs algorithm: Convergence and nonconvergence[END_REF] to cite just a few.

parameters ε ij , ∀i ∈ [1 . . . k], ∀j ∈ [1 . . . n]
One possibility is to perform weak-perspective iterations. Namely, the algorithm starts with a zero-distorsion weak-perspective approximation and then, at each iteration, it updates the perspective distorsions using eq. ( 33). To conclude, the robust perspective factorization algorithm can be summarized as follows:

Robust perspective factorization:

Initialization: Set ε ij = 0, ∀i ∈ [1 . . . k], ∀j ∈ [1 . . . n].
Use the same initialization step as the affine factorization with EM algorithm.

Iterate until convergence:

Affine factorization with EM: Iterate until convergence the E-and M-steps of the algorithm described in the previous section.

Euclidean upgrade: Recover the rotations, translations, and 3-D Euclidean coordinates from the affine shape and affine motion parameters.

Perspective update: Estimate new values for the parameters ε ij , ∀i ∈ Maximum a posteriori: After convergence choose in between inlier and outlier for each observation, i.e., max{α in ij ; α out ij }.

Multiple-camera calibration

In this section we describe how the solution obtained in the previous section is used within the context of multiple-camera calibration. As already described above, we are interested in the estimation of the external camera parameters, i.e., the alignment between a global reference frame (or the calibration frame) and the reference frame associated with each one of the cameras. We assume that the internal camera parameters were accurately estimated using available software.

There are many papers available that address the problem of internal camera calibration either from 3-D reference objects [START_REF] Faugeras | Three Dimensional Computer Vision: a Geometric Viewpoint[END_REF], 2-D planar objects [START_REF] Zhang | A flexible new technique for camera calibration[END_REF], 1-D objects [START_REF] Zhang | Camera calibration with one-dimensional objects[END_REF] or self-calibration, e.g., from point-correspondences [START_REF] Luong | Self-calibration of a moving camera from point correspondences and fundamental matrices[END_REF], [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], [START_REF] Luong | The Geometry of Multiple Images[END_REF].

Figure 2 shows a partial view of a multiple-camera setup as well as the one-dimensional object used for calibration. In practice we used three different camera configurations as depicted in Figure 4: two 30 camera configurations and one 10 camera configuration. These camera setups will be referred to as the Corner Case, the Arc Case, and the Semi-Spherical Case. Finding point correspondences accross the images provided by such a setup is an issue in its own right because one has to solve for a multiple wide-baseline point correspondence problem. We will briefly describe the practical solution that we retained and which maximizes the number of points that are matched over all the views. Nevertheless, in practice there are missing observations as well as badly detected image features, bad matches, etc. The problem of missing data has already been addressed. Here we concentrate on the detection and rejection of outliers.

We performed multiple camera calibration with two algorithms: The robust perspective factorization method previously described and bundle adjustment. We report a detailed comparison between these two methods. We further compare our robust method with a method based on M-estimators. As already mentioned, we use a simple 1-D object composed of four identical markers with known 1-D coordinates. These coordinates form a projective-invariant signature (the cross-ratio)

that is used to obtain 3-D to 2-D matches between the markers and their observed image locations.

With finely syncronized cameras it is possible to gather images of the object while the latter is freely moved in order to cover the 3-D space that is commonly viewed by all cameras. In the to strong perspective distorsion, that results in a number of connected components that cannot be easily matched with the four markers, and (c) only two connected components are visible in which case one cannot establish a reliable match with the four markers. In practice we perform a connected-component analysis that finds the number of blobs in each image. Each such blob is characterized by its center and second order moments, i.e., Figure 3 (d), (e), and (f). These blobs are matched with the object markers. In most of the cases the match is unambiguous, but in some cases a blob may be matched with several markers.

Let as before, s ij denote the center of a blob from image i that matches marker j. The second order moments of this blob can be used to compute an initial 2×2 covariance matrix C 0 ij for each such observation. Moreover, we introduce a binary variable, µ ij , which is equal to 0 if the observation s ij is missing and equal to 1 otherwise. The multiple-camera calibration algorithm This comparison requires the estimation of the rotations and translations allowing the alignment of the two reconstructed 3-D sets of points with the cameras. We estimate these rotations and translations using a set of control points. Indeed, both the robust perspective factorization and the bundle adjustment algorithms need a number of control points with known Euclidean 3-D coordinates. In practice, the calibration procedure provides such a set. This set of control points allows one to define a global reference frame. Let P c j denote the 3-D coordinates of the control points estimated with our algorithm, and let Q c j denote their 3-D coordinates provided in advance. Let λ, R, and t be the scale, rotation and translation allowing the alignment of the two sets of control points. We have:

min λ,R,t 8 j=1 λRQ c j + t -P c j 2 (36)
The minimizer of this error function can be found in closed form either with unit quaternions [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF] to represent the rotation R or with dual-number quaternions [START_REF] Walker | Estimating 3-d location parameters using dual number quaternions[END_REF] to represent the rigid motion R, t. Similarly, one can use the same procedure to estimate the scale λ ′ , rotation R ′ , and translation t ′ associated with the 3-D reconstruction obtained by bundle adjustment.

Finally, in order to evaluate the quality of the results we estimated the following measurements:

The 2D error is measured in pixels and corresponds to the RMS error between the observations and the predictions weighted by their posterior probabilities:

i,j α in ij s ij -ŝij (θ * ) 2 i,j α in ij 1/2
The 3D error is measured in milimeters and corresponds to the RMS error between the two sets of 3-D points obtained with our algorithm and with bundle adjustment:

n j=1 P j -P ′ j 2 n 1/2
The error in rotation is measured in degrees and depicts the average angular error of the rotation matrices over all the cameras. With the same notations as before, let R i and R ′ i be the rotations of camera i as obtained with our algorithm and with bundle adjustment. Let v be an arbitrary 3-D vector. The dot product (R i v) • (R ′ i v) is a reliable measure of the cosine of the angular discrepancy between the two estimations. Therefore the RMS error in rotation can be measured by the angle:

arccos   180 π v ⊤ k i=1 R ⊤ i R ′ i v k  
The error in translation is measured in milimeters with:

k i=1 t i -t ′ i 2 k 1/2
As already mentioned, we used three camera setups. All setups use identical 1024×768 Flea cameras from Point Grey Research Inc. 2 The intrinsic parameters were estimated in advance.

Two of the setups use 30 cameras, whereas the third one uses 10 cameras. We denote these setups as Corner Case, Arc Case and Semi-Spherical Case, based on the camera layout.

The results are summarized on Figures 4 and5 and on Tables 1 and2. Let us analyze in more detail these results. In the Corner Case there are 30 cameras and 292 3-D points. Hence, there are 8760 possible predictions out of which only 5527 are actually observed, i.e., 36% predictions correspond to missing 2-D data. The algorithm detected 5202 2-D inliers. An inlier is an observation with a posterior probability greater than 0.4. Next, the outliers are marked as missing data. Eventually, in this example, 285 3-D points were reconstructed (out of a total of 292) and all the cameras were correctly calibrated. The number of iterations of the robust perspective factorization algorithm (refered to as affine iterations) is equal to 7. On an average, there were 2.4 iterations of the EM algorithm. The obtained reconstruction has a smaller 2-D reprojection error (0.30 pixels) than the one obtained by bundle adjustment (0.58 pixels).

In any of the 3 calibration scenarios, the proposed method outperforms bundle adjustment results, as it can be observed in Table 2. This is in part due to the fact that the bundle adjustment algorithm does not have a mechanism for outlier rejection.

Figure 5 shows the evolution of the algorithm as it iterates from a weak-perspective solution to the final full-perspective solution. At convergence, the solution found by our method (shown in blue or dark in the absence of colors) is practically identical to the solution found by bundle adjustment (which is shown in grey).

Input Image 3-D Points Cameras

Corner Case Arc Case Semi-Spherical Case The bundle adjustment solution is shown in grey. The bundle adjustment solution does not perform outlier treatment.

3-D Reconstruction

The robust perspective factorization algorithm was also applied to the problem of 3-D reconstruction from multiple views. For this purpose we used images of objects using a single camera and a turning table. More specifically, we used the following data sets:

• The "Dino" and "Temple" data sets from the Middlebury's evaluation of Multi-View Stereo reconstruction algorithms;3

• The "Oxford dinausor" data set, 4 and

•
The "Square Box" data set.

We used the OpenCV5 pyramidal implementation of the Lucas & Kanade interest point detector and tracker [START_REF] Bouguet | Pyramidal implementation of the affine lucas kanade feature tracker -description of the algorithm[END_REF] to obtain the initial set of 2-D observations. This provides the 2k ×n measurement matrix S as well as the missing-data binary variables µ ij associated with each observation. note that these datasets are made available without the groundtruth 3-D data. They are typically used by the community to compare results for 3-D dense reconstructions. The rotation error stays within 3 degrees. The translation error is very small because, in this case we aligned the camera centers and not the 3-D coordinates of some reconstructed points.

The results obtained with the Oxford "Dinausor" need some special comments. Because of the very large percentage of missing data, we have been unable to initialize the solution with the PowerFactorization method. Therefore, we provided the camera calibration parameters for initialization. However, this kind of problem can be overcome by using an alternative affine factorization algorithm [START_REF] Ph | Algorithms for batch matrix factorization with application to structure-from-motion[END_REF].

In order to further assess the quality of our results, we used the 3-D reconstructed points to build a rough 3-D mesh and to further apply a surface-evolution algorithm to the latter in order to obtain a more accurate mesh-based 3-D reconstruction [START_REF] Zaharescu | Transformesh: a topology-adaptive mesh-based approach to surface evolution[END_REF]. The results are shown on Figure 8. 6

Comparison with other methods

As already mentioned in section 3, our robust ML estimator has strong similarities with Mestimators and their practical implementation, i.e., IRLS [START_REF] Stewart | Robust parameter estimation in computer vision[END_REF]. Previous work on robust affine factorization has successfully used the following reweighting function φ that corresponds to the truncated quadratic:

φ(x) =    1 if |x| < k k 2 x 2 otherwise (37)
It is therefore tempting to replace the EM procedure of our algorithm with an IRLS procedure, which amounts to replace the posterior probabilities of inliers α in ij given by eq. ( 10) with the weights φ ij given by eq. [START_REF] Sturm | A factorization based algorithm for multi-image projective structure and motion[END_REF]. The latter tends to zero most quickly allowing aggressive rejection of outliers. One caveat is that the efficiency of IRLS depends on the tuning parameter k.

Unfortunately the latter cannot be estimated within the minimization process as is the case with the covariance matrix. However, we noted that the results that we obtained do not depend on the Dino Temple Box Dinausor Figure 8: Dense reconstruction results: A rough mesh obtained from the 3-D reconstructed points (middle) and the final dense reconstruction (bottom) after surface evolution using the method described in [START_REF] Zaharescu | Transformesh: a topology-adaptive mesh-based approach to surface evolution[END_REF].

choice of k. In all the experiments reported below, we used k = 1, The plot of the truncated quadratic for different k values is plotted in Figure 9. 4.

In the Corner case the quality of the results are very similar: our algorithm accepted 94% of the total number of observations as inliers and reconstructed 99.3% of the total number of 3-D points, while IRLS accepted all the observations as inliers and reconstructed all the 3-D points.

Similar results are obtained in the Arc and Semi-Spherical cases, where the proposed method performs slightly better. Both algorithms were able to reconstruct the Dino and the Temple, but our algorithm yields more accurate results. Outlier detection is summarized in Table 5.

A more thorough comparison with robust as well as non robust 3-D reconstruction methods is provided in Figure 10. The proposed algorithm is denoted by "Persp. Power Factorization (Bayesian)", while the IRLS method is named "Persp. Power Factorization (IRLS -Truncated Quadratic)" and the non-robust method is called "Persp. Power Factorization (Not Robust)". Affine factorization algorithms are also presented, together with the results of bundle adjustment. The bundle adjustment method was always initialized using the PowerFactorization method. The robust perspective factorization method proposed in this paper is the most resilient

Conclusions

In this paper we described a robust factorization method based on data clustering and on the EM algorithm. First we recalled the classical maximum-likelihood approach within which all the observations are supposed to be independent and identically distributed. This amounts to classify We performed extensive experiments with two types of data sets: multiple-camera calibration and 3-D reconstruction. We compared the calibration results of our algorithm with the results obtained using other methods such as the bundle adjustment technique and IRLS. It is interesting to notice that there is almost no noticeable quantitative difference between our algorithm and a nonlinear optimization method such as bundle adjustment. The 3-D reconstruction results obtained with a single camera and objects lying on a turntable are also very good. Whenever possible, we compared our results with ground-truth data, such as the external camera parameters provided by the Middlebury multi-view stereo data set. In order to further assess the 3-D reconstruction results, we used the output of the robust perspective factorization method, namely a cloud of 3-D points, as input of a mesh-based reconstruction technique.

Our Gaussian/uniform mixture model and its associated EM algorithm may well be viewed as a robust regression method in the spirit of M-estimators. We compared our method with IRLS using a truncated quadratic loss function. The results show that our method performs slightly better, although we believe that these results are only preliminary. A thorough comparison between outlier detection using probability distribution mixture models on one side, and robust loss functions on the other side is a topic in its own right. In the future we plan to extend our method to deal the more difficult problem of multiple-body factorization.

A Derivation of equation (19)

We recall eq. ( 17):

Q M L = 1 2 i,j (s ij -ŝij (θ)) ⊤ C -1 (s ij -ŝij (θ)) + log(det C)
Taking the derivative with respect to the entries of the 2×2 matrix C we obtain:

∂Q M L ∂C = 1 2 i,j -C -⊤ (s ij -ŝij (θ))(s ij -ŝij (θ)) ⊤ C -⊤ + C -⊤ (38) = -C -⊤ 1 2 i,j (s ij -ŝij (θ))(s ij -ŝij (θ)) ⊤ C -⊤ + m 2 C -⊤ (39)
where m = k × n. By setting the derivative to zero we obtain eq. ( 19):

C = 1 m i,j (s ij -ŝij (θ))(s ij -ŝij (θ)) ⊤

  are available, then this corresponds to a weakperspective camera model that is closer to the true perspective model. If the true values of the perspective-distortion parameters are available, the corresponding weak-perspective model corresponds exactly to the perspective model. Hence, the problem reduces to affine factorization followed by Euclidean upgrade. Numerous iterative algorithms have been suggested in the literature for estimating the perspective-distortion parameters associated with each 2-D observation,

[ 1 .

 1 . . k], ∀j ∈ [1 . . . n]. If the current depth values are identical with the previously estimated ones, then terminate, else iterate.

Figure 2 :

 2 Figure 2: (a): Partial view of a 30-camera setup. (b): The calibration data are gathered by moving a one-dimensional object in the common field of view of the cameras.

  three examples below we used 73, 58, and 32 frames, i.e., 292, 232, and 128 3-D points. The number of cameras in each setup is 30, 30 and 10, respectively. Therefore, in theory there should be 8760, 6960, and 1280 2-D observations.

Figure 3 :

 3 Figure 3: Top: These are typical images where the number of connected components depend on the position and orientation of the calibrating object with respect to the cameras. Bottom: Detected blobs with their centers and associated covariance, i.e., second-order moments.

Figure 3

 3 Figure 3 depicts three possible image configurations: (a) four distinct connected components that correspond without ambiguity to the four markers, (b) a degenerate view of the markers, due

Figure 4 :

 4 Figure 4: Multiple camera calibration results. Left: A typical set of 2-D observations associated with one camera. Middle: Reconstructed 3-D points with our method (blue) and with bundle adjustment (grey). Right: Camera calibration results obtained with our method (blue) and with bundle adjustment (grey) .

Figure 5 :

 5 Figure 5: Iterations of the robust perspective factorization algorithm in the Arc Case and comparison with bundle adjustment. The first iteration corresponds to weak-perspective factorization.

Fig- ures 6 Figure 6 :Figure 7 :

 667 Figure 6: Groundtruth data is represented in gray (light) colour, whereas reconstruction results are represented in blue (dark) colour.

Figure 9 :

 9 Figure 9: Reweighting function for k = 0.1, 0.5, 1, 2, 5, 10.

Figure 11

 11 Figure 11 compares our method (a), with the bundle adjustment method (b), in the Arc Caseand when 20% of the input data was corrupted by high-amplitude noise (σ = 0.20 of the image size). On both figures the ground truth is shown in grey and the result of the algorithm is shown in blue (or dark in the absence of color). Notice that with this level of data perturbation,

Figure 10 :

 10 Figure 10: Behavior of various robust and non robust algorithms when an increasing percentage of the input data are corrupted by high-amplitude noise, namely σ = 0.2 of the image size.

  (a) Robust perspective factorization (b) Bundle Adjustment

Figure 11 :

 11 Figure 11: Calibration results in the Arc Case for (a) the proposed method and for (b) bundle adjustment method, when 20% of the input data are corrupted with high-amplitude noise, namely σ = 0.2 of the image size. The 2-D reprojection error is of 0.71 pixels for (a) and 15.84 pixels for (b). The groundtruth is represented in gray.

Table 1 :

 1 ) . Summary of the camera calibration results for the three setups.

	Multi-camera calibration	Corner Case Arc Case Semi-Spherical Case
	Input # Cameras		30	30	10
	# 3-D Points		292	232	128
	# 2-D Predictions		8760	6960	1280
	# Missing observations		36%	0%	33%
	# 2-D Observations		5527	6960	863
	Results # 2-D Inliers		5202	6790	784
	# 3-D Inliers		285	232	122
	2D error (pixels)		0.30	0.19	0.48
	3D error (mm)		6.91	2.65	4.57
	Rot. error (degrees)		0.13	0.18	0.27
	Tr. error (mm)		27.02	9.37	24.21
	# Aff. iter. (# EM iter.)		7 (2.4)	11 (2)	8 (3.2)
	Multi-camera calibration Corner Case Arc Case Semi-Spherical Case
	Proposed Method -2D error (pixels)	0.30	0.19	0.48
	Bundle Adjustment -2D error (pixels)	0.58	0.61	0.95

Table 2 :

 2 Comparison between the proposed method and bundle adjustment for the three camera calibration setups.

	1 st iteration	2 nd iteration	3 rd iteration

2D Err. = 59.32 2D Err.= 18.43 2D Err. =15.92 4 th iteration 5 th iteration 8 th iteration 2D Err. = 4.53 2D Err. = 1.61 2D Err. = 0.24

Table 3 :

 3 Summary of the 3-D reconstruction results for the four data sets.

	3-D reconstruction	Dino	Temple	Box	Dinausor
	Input # Views	12	16	64	36
	Size of S matrix	24×480 32×758 128×560 72×1516
	# 2-D predictions	5760	12128	35840	54576
	% Missing observations	11%	21%	17%	85%
	# 2-D Observations	5140	9573	29733	8331
	Results # 2-D Inliers	3124	6811	25225	7645
	# 3-D Inliers	370	720	542	1437
	2D error (pixels)	0.82	0.93	0.69	0.33
	Rot. error (degrees)	1.49	2.32	-	0.00
	Trans. error (mm)	0.01	0.01	-	0.00
	# Aff. iter. (# EM iter.)	7 (4)	7 (3.14)	9 (3)	7 (3.29)

Table 4 :

 4 Dataset Method 2-D Inliers 3-D Inliers 2-D err. 3-D err. Rot. err. Trans. err. Comparison between robust perspective factorization results using EM and IRLS. The figures in paranthesis correspond to the total number of observations (third column) and to the total number of expected 3-D points (fourth column).

	Corner	EM	5202 (5527) 285 (292)	0.30	6.91	0.13	27.02
		IRLS 5526 (5527) 288 (292)	0.40	6.91	0.14	26.61
	Arc	EM	6790 (6960) 232 (232)	0.19	2.65	0.18	9.37
		IRLS 6960 (6960) 232 (232)	0.22	2.54	0.16	8.78
	Semi-Spherical	EM	784 (863)	122 (128)	0.48	4.57	0.27	24.21
		IRLS	862 (863)	128 (128)	0.62	4.66	0.29	23.91
	Dino	EM	3124 (5140) 370 (480)	0.82	-	1.49	0.01
		IRLS 3411 (5140) 390 (480)	2.57	-	2.13	0.01
	Temple	EM	6811 (9573) 720 (758)	0.93	-	2.32	0.01
		IRLS 7795 (9573) 731 (758)	1.69	-	2.76	0.03
			Corner Arc Semi-Spherical Dino Temple
		EM	6%	2%	9%	39%	29%
		IRLS	0%	0%	0%	34%	19%

Table 5 :

 5 Percentage of outliers detected by the two algorithms.to high-amplitude noise. It generally performs better than the IRLS method and provides a clear advantage against the non-robust methods, which exit the graphs as soon as the noise level increases. As it can be observed, in the Semi-Spherical Case, the solution deteriorates a lot faster in the presence of noise, due to the lack of the redundancy in the data (128 3-D points and 10 cameras, versus 292 points and 30 cameras in the Corner Case and 232 points and 30 cameras in the Arc Case).

P (z ij = inlier|s ij )P (s ij ) = P (s ij |z ij = inlier)P (z ij = inlier)

http://www.ptgrey.com/products/flea/

http://vision.middlebury.edu/mview/data/

http://www.robots.ox.ac.uk/~vgg/data/data-mview.html

http://www.intel.com/technology/computing/opencv/

They are also available at http://vision.middlebury.edu/mview/eval/.

can now be paraphrased as follows:

Multiple camera calibration:

Initialization: Use eq. ( 18) to estimate the affine shape and motion parameters in the presence of some missing data:

Estimate the initial covariance matrix using eq. ( 19):

Iterate until convergence:

Affine factorization with EM: Iterate until convergence the E-and M-steps of the algorithm described in section 3.

Euclidean upgrade: Recover the rotations, translations, and 3-D Euclidean coordinates from the affine shape and affine motion parameters. Both the initialization and the M steps of the above algorithm perform affine factorization in the presence of uncertainty and missing data. In [START_REF] Zaharescu | Multiple camera calibration using robust perspective factorization[END_REF] we compared several such algorithms and we came to the conclusion that the PowerFactorization algorithm outperforms the other tested algorithms. In order to assess quantitatively the performance of our algorithm, we compared it with an implementation of the bundle adjustment method along the lines described in [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF].

B Derivation of equation (21)

When considering isotropic covariance, C = σ 2 I, hence det C = σ 4 , and the equation becomes:

By taking the derivative with respect to σ 2 , we obtain:

By setting the derivative to zero, we obtain eq. ( 21):

where m = k × n.