Adrien Bartoli

A Unified Framework for Quasi-Linear Bundle Adjustment

Obtaining 3D models from long image sequences is a major issue in computer vision. One of the main tools used to obtain accurate structure and motion estimates is bundle adjustment. Bundle adjustment is usually performed using nonlinear Newton-type optimizers such as Levenberg-Marquardt which might be quite slow when handling a large number of points or views.

We investigate an algorithm for bundle adjustment based on quasi-linear optimization. The method is straightforward to implement and relies on solving weighted linear systems obtained as simple functions of the input data. Important features are that (i) the original cost function is preserved, (ii) different projection models, either calibrated or not, are handled in a unified framework and (iii) any number of views and points as well as missing data can be handled.

Experimental results on simulated and real data show that the algorithm is as accurate as standard techniques while requiring less computational time to converge.

Introduction

Recovering metric structure and motion from uncalibrated images is a central problem for computer vision. Most of the time, a sub-optimal solution is obtained in projective space for motion, then for structure [START_REF] Beardsley | 3D model acquisition from extended image sequences[END_REF], or jointly [START_REF] Sturm | A factorization based algorithm for multi-image projective structure and motion[END_REF]. Selfcalibration is then performed to upgrade the reconstruction to metric space.

In this paper, we investigate a very simple technique for bundle adjustment. Experimental results reveal that it is as accurate as traditional techniques while requiring less CPU time. We consider both the full perspective and the affine projection, either calibrated or uncalibrated. We may end up therefore with projective, affine or metric reconstructions.

Let us consider point features. Bundle adjustment basically consists in finding the minimum of a cost function defined as the sum of discrepancies between actual and predicted image points. These discrepancies are measured using the Euclidean distance. The optimization is performed over a set of parameters that represents structure and motion using most of the time non-linear Newton-type optimization techniques such as Levenberg-Marquardt [START_REF] Triggs | Bundle ajustment -a modern synthesis[END_REF]. This technique has the following drawbacks. Firstly, handling the free scale of homogeneous entities such as vectors or matrices representing points or cameras may be non-trivial. Secondly, it requires the computation of at least the first order Jacobian matrix of the residuals with respect to structure and motion parameters which might be non-trivial. Thirdly, the computational cost may be non-negligible since the Hessian matrix has to be inverted, even if specific techniques have been proposed to speed up the process [START_REF] Triggs | Bundle ajustment -a modern synthesis[END_REF].

In this paper, we adapt the resection-intersection technique [START_REF] Chen | Efficient iterative solution to M-view projective reconstruction[END_REF][START_REF] Mahamud | Provablyconvergent iterative methods for projective structure and motion[END_REF][START_REF] Triggs | Bundle ajustment -a modern synthesis[END_REF] to perform bundle adjustment using quasilinear optimizations. Our contribution differs from previous ones in that (i) we use the original cost function of bundle adjustment, which preserves optimality and (ii) we handle a great variety of camera models in a unified manner. In particular, we deal with calibrated configurations where nonlinear constraints hold on the recovered motion. The final algorithm is simplicity itself. The method relies on rewriting the non-linear Euclidean distance used in bundle adjustment as a weighted bilinear algebraic distance [START_REF] Hartley | Minimizing algebraic error[END_REF]. In practice, it consists in iteratively solving weighted least squares systems that are simple functions of image point positions. Experimental results reveal that it performs as well as standard techniques in terms of convergence accuracy while greatly reducing computational cost.

This paper is organized as follows. In §2, we give our notation and preliminaries. §3 introduces the principle and the algorithm for quasi-linear bundle adjustment. These results are independent of the camera model. The approach is specialized in § §4.1 and 4.2 to uncalibrated and calibrated perspective cameras respectively. Finally, §5 validates the approach on simulated and real data and §6 gives our conclusions and perspectives.

Notation and preliminaries.

Vectors are typeset using bold letters (q, Q), matrices using sans-serif fonts (P, S) and scalars in italics (x, w). We make no formal distinction between coordinate vectors and physical entities. Everything is expressed in homogeneous coordinates, e.g. q T ∼ (q 1 q 2 q 3 ) is an image point, where ∼ means equality up to scale and T is vector/matrix transpo-sition.

The standard cost function for bundle adjustment is defined by C = i j d 2 (q ij , q ij ), where d 2 is the squared Euclidean distance, q ij the reprojection of the i-th estimated 3D point in the j-th image and q ij the corresponding image point. We drop the indices i and j for clarity. Bundle adjustment can be formulated as min P,Q C, where P and Q respectively designate the estimates of the parameters of motion and structure.

Distances in the image. The Euclidean distance is employed for bundle adjustment and in many criteria since it is physically meaninful. Most of the time, its non-linearity induces that of the corresponding criterion which therefore can not be used to linearly estimate e.g. structure or motion. On the other hand, consider the squared algebraic distance d a defined by [START_REF] Hartley | Minimizing algebraic error[END_REF]:

d 2 (x, y) = w 2
x,y d 2 a (x, y) where w x,y = (x 3 y 3 ) -2 . (1) The non-linearity is hidden in the weight factors w x,y which can be thought of as a bias of the algebraic distance with respect to the Euclidean one. A convenient expression may be derived as follows:

d 2 a = ||S[x] × y|| 2 and S = -1 0 0 0 1 0 , (2)
where ||.|| 2 is the L 2 or squared Frobenius norm for vectors or matrices respectively and [.] × the skew-symmetric 3 × 3matrix associated with the cross product, i.e.

[v] × q = v × q. The algebraic distance does not have any physical meaning since it depends upon the relative scale of the homogeneous representation of x and y. Its advantage is that it is bilinear. Since in most estimation problems one of the points compared is known, d a becomes linear for the other.

Quasi-linear optimization. The relationship (1) between the Euclidean and the algebraic distance is the basis for quasilinear optimization. The principle is to initialize weight factors to unity and compute a biased estimate using the weighted algebraic distance. When such an estimate has been obtained, weight factors can be computed using [START_REF] Beardsley | 3D model acquisition from extended image sequences[END_REF]. The process is then iterated until convergence. Typically, 3 or 4 iterations are enough.

Quasi-Linear Bundle Adjustment

We derive the bundle adjustment algorithm, which is valid for most camera models. We then detail camera-dependent elements.

Consider the initial bundle adjustment problem which involves optimizing simultaneaously over the motion and structure parameters. It can be split, as indicated in [START_REF] Chen | Efficient iterative solution to M-view projective reconstruction[END_REF][START_REF] Mahamud | Provablyconvergent iterative methods for projective structure and motion[END_REF][START_REF] Zhang | Motion and structure from two perspective views: From essential parameters to euclidean motion via fundamental matrix[END_REF] and becomes min P (min Q C). This is strictly equivalent to the initial problem in that the same minima are kept. We split accordingly the optimization process into two distinct steps, for motion on the one hand and for structure on the other hand.

These optimization steps can be alternatively performed to refine structure and motion. This scheme is termed resectionintersection. Each step can then be conducted independently using non-linear optimizers such as Levenberg-Marquardt, as proposed in [START_REF] Zhang | Motion and structure from two perspective views: From essential parameters to euclidean motion via fundamental matrix[END_REF]. The advantages of doing so are numerous:

• it reduces the complexity, since the optimization consists of (m -1) independent outer steps1 over an 11dimentional space and n independent 3-dimensional inner steps (m is the number of views of n is the number of points); • it performs as well as directly optimizing over all the parameters in terms of convergence accuracy, as reported by [START_REF] Chen | Efficient iterative solution to M-view projective reconstruction[END_REF][START_REF] Mahamud | Provablyconvergent iterative methods for projective structure and motion[END_REF][START_REF] Zhang | Motion and structure from two perspective views: From essential parameters to euclidean motion via fundamental matrix[END_REF] and confirmed by our results; • each of the two alternated optimization steps can be performed quasi-linearly, see the next section.

The algorithm, given below, needs an initialization of either structure or motion whereas standard bundle adjusters require both:

1. initialize weight factors to unity: w q, q = 1. These are not necessary for the affine camera model; 2. quasi-linearly solve for structure; 3. quasi-linearly solve for motion; 4. iterate steps 2 and 3 until convergence (see below);

The estimation of weight factors, which corrects the bias of the algebraic distance with respect to the Euclidean one, is done inside steps 2 and 3.

Convergence.

Convergence is achieved when steps 2 and 3 become idempotent for structure and motion, i.e. when the difference between two consecutive values of the cost function C is lower than a threshold, typically 10e -8 pixel, or when the residual error begins to increase.

Resection-Intersection

Uncalibrated Pin-Hole Cameras

We consider uncalibrated full perspective projection to model cameras. This yields projective reconstructions. We derive the equations that lead to structure, then motion, quasilinear estimation. Consider one term d 2 (q, q) of the double sum of the cost function C. Using the link (1) between the Euclidean and the algebraic distances d and d a , then the expression (2) of the algebraic distance, each of these terms can be written d 2 (q, q) = w 2 q, q ||S[q] × q|| 2 . Using perspective projection q ∼ PQ and some minor algebraic manipulations:

d 2 (q, q) = w 2 q, q p T 2 Q -q 2 p T 3 Q p T 1 Q -q 1 p T 3 Q 2 , (3)
where p T i are the rows of P. Using (1), we obtain the expression for the weight factors as:

w q, q = (p T 3 Q) -2 . ( 4
)
Intersection -Solving for the structure. Using equation (3), and if weight factors w q, q are assumed to be known, one can recover the structure by solving the linear system M us Q = 0 for each point feature Q considered where:

(M us ) 2m×4 =        . . . w q q    p T 2 -q 2 p T 3 p T 1 -q 1 p T 3    . . .        . (5)
Each of the m views where the point Q is visible provides 2 equations. The solution is the unit singular vector corresponding to the smallest singular value of M us which can be obtain using e.g. SVD. This naturally ensures the free scale of the recovered homogeneous vector since ||Q|| 2 = 1.

The process is then iterated while re-estimating weight factors w q, q to reduce the bias of the algebraic distance until convergence, which is assessed as in §3.

Resection -Solving for the motion. We use equation ( 3) to form a linear system M um p that can be solved using SVD where p ∼ (p T 1 p T 2 p T 3 ) encapsulates the rows of P and:

(M um ) 2n×12 =        . . . w q q    0 T 4 Q T -q 2 Q T Q T 0 T 4 -q 1 Q T    . . .       
.

Each of the n points that project onto view P provides 2 equations. The iteration is then conducted while re-estimating the weight factors, as for the structure. The free scale of P is fixed since

||P|| 2 = ||p|| 2 = 1.

Calibrated Pin-Hole Cameras

We consider the calibrated full perspective camera model. This leads to metric reconstructions. In this case, resection is termed pose estimation and has been widely studied, see e.g. [START_REF] Haralick | Pose estimation from corresponding point data[END_REF]. When the motion is frozen, solving for the structure can be conducted using the quasi-linear algorithm of §4.1.

Pose estimation from minimal data, i.e. n = 3 points is non-linear by nature since there exist four solutions in this case. However, when an initial guess is known, it can be refined using linear algebra. Difficulties arise because of nonlinear constraints on the orthonormality of the recovered rotation matrix. An algorithm for pose refinement with n ≥ 3 points that fits into our quasi-linear framework is described in the next section.

Quasi-linear pose estimation. Let P ∼ K(R t) be the projection matrix where K is a 3×3 matrix containing the known intrinsic parameters, R a 3 × 3 rotation matrix and t a 3vector. The unknowns are therefore R and t. Under the assumption of zero skew and unity aspect ratiowhich are valid for most modern cameras, distances measured onto the image plane are equal to distances measured onto the retina up to a scale factor as d 2 (q, q) = f 2 d 2 (x, x) where x = K -1 q and x -1 = K -1 q are retinal points and f is the focal length (in number of pixels).

Linearly update the motion. Each term d 2 (q, q) of the double sum of the cost function C can be expanded as in the uncalibrated case, but using retinal points. This gives

d 2 (q, q) = f 2 w 2 x, x ||S[x] × x|| 2 . Using then perspective pro- jection x ∼ (R t)Q of the point Q T ∼ ( QT α) on the retina, we obtain d 2 (q, q) = f 2 w 2 x, x ||S[x] × (R Q + αt)|| 2 .
If we approximate the rotation between two iterations of the algorithm by a first order taylor expansion, (I + [w] × ), we obtain after some minor algebraic manipulations d

2 (q, q) = f 2 w 2 x, x ||S[x] × R Q -S[x] × [ Q] × w + αS[x] × t|| 2 + O,
where O contains higher order terms. By neglecting O, we can estimate the first order pose update m T = (w T t T ) by solving the linear system M cm m = d cm , where:

(M cm ) 2n×6 =    . . . f 2 w 2 x, x S[x] × R[ Q] × -αS[x] × . . .    (d cm ) 2n =    . . . f 2 w 2 x, x S[x] × R Q . . .    .
Algorithm. The algorithm for quasi-linear pose estimation combines first order updates of pose and re-estimation of calibrated weight factors w x, x = 1 f w q, q : 1. solve for pose update using M cm m = d cm ; 2. update rotation: R ← RR w , where R w gives the rotation corresponding to w T = (w 1 w 2 w 3 ) using e.g. Euler angles; 3. re-estimate calibrated weight factors w x, x using (4); 4. iterate steps 1-3 until convergence, see §3.

Experimental Results

In the following two sections, we compare previously described algorithms to existing ones. We use simulated data to compare them in the most general case, i.e. projective. We then consider real images and both projective and metric spaces. More detailled results, in particular on pose estimation and for the affine camera will be made available in a longer version of the paper. We compare the following algorithms. LIN is structure and motion initialization [START_REF] Beardsley | 3D model acquisition from extended image sequences[END_REF], QLIN is the proposed quasi-linear method and LM relies on the Levenberg-Marquardt non-linear optimizer. Note that our implementation uses analytical differentiation and a sparse inverser for the Hessian matrix, see [START_REF] Triggs | Bundle ajustment -a modern synthesis[END_REF]. We also compared a full Hessian inversion-based implementation that we denote as LM-full.

The threshold for convergence (see §3) is the same for QLIN and LM and is equal to 10e -8 . Image point coordinates are standardized such that they lie in [-1 . . . 1] 2 .

Simulated Data

The simulated scene consists of 3D point features chosen at random inside a cube of 2 meters side length. We simulate two camera configurations. The first one is stable and the second one contains nearly unreconstructable points. We measure the reprojection error to assess the convergence accuracy and the CPU time to convergence to compare computational costs. In order to simulate realistic situations, we adopted the following parameters: the focal length of the cameras is 1000 (expressed is number of pixels) and the image size is 512 × 512. Values plotted have been averaged over 50 trials. Median values gave graphs similar to those for the means. The curves for 3D errors show expected results, see e.g. [START_REF] Sturm | A factorization based algorithm for multi-image projective structure and motion[END_REF].

Convergence accuracy. The standard configuration is 50 points, 10 views and an added image noise with a Gaussian distribution of 0.5 pixel standard deviation.

Figures 1 (a) shows the results when varying the image noise. We observe that for camera configuration 1, there is no difference between QLIN and LM, even when the starting point given by LIN lies quite far from the final result. The error for these two methods degrades gracefully as image noise increases. For camera configuration 2, figure 1 (a) shows that there is a slight difference between QLIN and LM only when the level of added image noise is beyond 1 pixel. LM-full, not shown on these graphs, gave exactly the same results as LM. Computation time. We observe on figure 1 (b) that the quasi-linear optimization requires less CPU time to convergence than the non-linear one whatever the number of views.

-conf. 1 QLIN -conf. 1 LM -conf. 1 LIN -conf. 2 QLIN -conf. 2 LM -conf. 2 (a) (b)
Results from another experiment where added image noise is varied from 0 to 2 pixels, shows that CPU time to convergence linearly increases with the added image noise. The CPU time needed by LM is twice thus needed by QLIN. LM-full gave results out of range of these graphs.

Real Images

We use the 181 images of the hotel sequence 2 and 197 corner correspondences. These corners have been automatically detected and semi-automatically matched through the different views. Calibration data were approximately known for this sequence. The results given in table 1 show that QLIN converges to the same residual error as LM and LM-full. In terms of computational cost, we observe that QLIN is roughly 2.5 times faster than LM.

Conclusions and Perspectives

We addressed the problem of structure and motion refinement using bundle adjustment. We perform this task based on quasi-linear optimization while keeping the original cost function of bundle adjustment. We split the problem into structure refinement on the one hand and motion refinement on the other hand. The principle is then to rewrite the Euclidean distance used in bundle adjustment as a weighted algebraic distance. Quasi-linear optimization iteratively updates these weights while refining structure or motion.

The result is a disarmigly simple algorithm that comes out to be a loop over two weighted linear systems constructed as simple functions as the input data.

We conducted numerous experiments on simulated and real data. In the light of these results, we observed that this algorithm is as accurate as standard Levenberg-Marquardtbased bundle adjusters in terms of convergence accuracy while being much faster in terms of computational cost.

Figure 1 .

 1 Figure 1. (a): reprojection error at convergence versus added image noise. The curves for QLIN and LM and the first configuration, and for LM and the second configuration are nearly undistinguishable. (b): CPU time to convergence versus the number of views.

2

  These data have been provided by the Modeling by Videotaping group in the Robotics Institute, Carnegie Mellon University.

	step	approach	first 5 images rep. error (pixel) CPU time (second)	full 181 images sequence rep. error (pixel) CPU time (second)
	init.	-	9.4209	-	6.2524	-
		LM-full	0.0357	1056 ≈17 mn.	-	-
	ba	LM	0.0357	2.55	2.2026	96
		QLIN	0.0391	0.95	2.2034	37

Table 1 . Reprojection errors (pixel) and CPU time to con- vergence (second) for various algorithms and sequence lengths.

 1

we do not optimize the first camera parameters to (partially) fix the coordinate frame.

Acknowledgement. I wish to thank Frederik Schaffalitzky for providing his implementation of LM used in §5.


Adrien Bartoli

On the Non-Linear Optimization of Projective Motion Using Minimal Parameters

I address the problem of optimizing projective motion over a minimal set of parameters. Most of the existing works overparameterize the problem. While this can simplify the estimation process and may ensure well-conditioning of the parameters, this also increases the computational cost since more unknowns than necessary are involved. I propose a method whose key feature is that the number of parameters employed is minimal. The method requires singular value decomposition and minor algebraic manipulations and is therefore straightforward to implement. It can be plugged into most of the optimization algorithms such as Levenberg-Marquardt as well as the corresponding sparse versions. The method relies on the orthonormal camera motion representation that I introduce here. This representation can be locally updated using minimal parameters. I give a detailled description for the implementation of the two-view case within a bundle adjustment framework, which corresponds to the maximum likelihood estimation of the fundamental matrix and scene structure. Extending the algorithm to the multiple-view case is straightforward. Experimental results using simulated and real data show that algorithms based on minimal parameters perform better than the others in terms of the computational cost, i.e. their convergence is faster, while achieving comparable results in terms of convergence to a local optimum. An implementation of the method will be made available.

Introduction

The problem of recovering structure and motion from images is one of the central challenges for computer vision. The use of image feature correspondences (e.g. points, lines) through different views and the study of geometrical aspects of the image formation process have led to numerous techniques acting in metric, affine or projective space, depending on whether camera calibration is fully, partially or not available.

Most of the time, a sub-optimal solution is obtained using linear techniques for motion then for structure recovery [START_REF] Beardsley | 3D model acquisition from extended image sequences[END_REF] or jointly [START_REF] Sturm | A factorization based algorithm for multi-image projective structure and motion[END_REF] and subsequently refined. One of the most efficient techniques for such a structure and motion optimization is bundle adjustment. It involves minimizing a non-linear cost function based on the discrepancy between reprojected and original image features [START_REF] Slama | Manual of Photogrammetry, Fourth Edition. American Society of Photogrammetry and Remote Sensing[END_REF][START_REF] Triggs | Bundle ajustment -a modern synthesis[END_REF]. The behaviour of such techniques, in terms of convergence properties to a local optimum of the cost function and computational cost, greatly depends on the algebraic representation of the problem, i.e. of structure and motion, and in particular, numerical conditioning and whether the number of parameters employed is minimal. For bundle adjustment, preserving the original noise model on image features is also crucial.

While bundle adjustment is theoretically well-defined (it corresponds to the maximum likelihood estimator, see e.g. [START_REF] Triggs | Bundle ajustment -a modern synthesis[END_REF]) there are still practical optimization problems since the employed cost functions have many local minima [START_REF] Oliensis | The error surface for structure and motion[END_REF] where optimization processes may get trapped.

I address the problem of representing motion. The goal is to obtain minimal estimators, i.e. where the number of parameters considered for optimization is minimal. I focus on the projective case, i.e. when camera calibration is unknown. This topic is of primary importance since an accurate projective reconstruction is necessary to subsequently succeed in self-calibration.

Consider two perspective views of a rigid scene. The image signature of their relative position is the projective two-view motion or the epipolar geometry, usually described by the (3 × 3) rank-2 and homogeneous fundamental matrix [START_REF] Luong | The fundamental matrix: Theory, algorithms and stability analysis[END_REF][START_REF] Zhang | Determining the epipolar geometry and its uncertainty: A review[END_REF]. A fundamental matrix has 7 degrees of freedom. Therefore, 7 parameters should be enough to optimize the projective two-view motion.

However, it has been seen that there does not exist a universal parameterization of the fundamental matrix using 7 parameters. This is due to the non-linear rank-2 constraint and the free-scale ambiguity. Existing works may fall into the following categories:

-overparameterization, e.g. the 12 entries of a perspective camera matrix [START_REF] Hartley | Projective reconstruction and invariants from multiple images[END_REF]. More unknowns than necessary are involved to simplify the representation; -multiple minimal parameterizations, 3 in [START_REF] Bartoli | Projective structure and motion from two views of a piecewise planar scene[END_REF] or 36 in [START_REF] Zhang | Determining the epipolar geometry and its uncertainty: A review[END_REF]; -minimal parameterizations combined to image transformations [START_REF] Bartoli | Projective structure and motion from two views of a piecewise planar scene[END_REF][START_REF] Zhang | Estimating the fundamental matrix by transforming image points in projective space[END_REF] to reduce the number of parameterizations.

Other techniques optimize over the 9 parameters of the fundamental matrix while adding the non-linear rank-2 constraint and the normalization constraint as virtual measurements. This paper makes the following contributions. Firstly, I address the projective two-view motion case in §2. I introduce what I call the orthonormal representation of projective two-view motion. Based on this, I show how one can non-linearly estimate the projective two-view motion using a minimal number of 7 parameters. An important point is that this method does not depend upon the optimization technique considered.

Secondly, I illustrate the use of this method in a bundle adjustment framework based on [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] in §3. The result is a minimal maximum likelihood estimator for the fundamental matrix as well as scene structure. The reader who is interested into practical issues only should refer directly to this section.

Thirdly, I extend the framework to multiple views in §4 where I introduce the orthonormal representation of projective multiple-view motion. I derive, similarly to the two-view case, a means to perform optimization over a minimal set of parameters.

Finally, experimental results on simulated and real data are shown in § §5 and 6 respectively. They show that algorithms based on minimal motion parameters, and in particular on the orthonormal representation, perform better than the others, in terms of computational cost, while achieving equivalent performances in terms of convergence properties. This is followed by my conclusions in §7.

The Projective Two-View Motion

Preliminaries

Let us consider two (3 × 4) uncalibrated perspective camera matrices. Due to homogeneity, each one has 11 degrees of freedom. Since there is a 15-degrees of freedom coordinate frame ambiguity on structure and motion, the projective two-view motion has 11 • 2 -15 = 7 degrees of freedom.

The corresponding (3 × 3) fundamental matrix F has 9 entries but 7 degrees of freedom since it is homogeneous and has rank 2. It allows one to extract projection matrices for the two views while fixing the coordinate frame. These projection matrices constitute a realization of the fundamental matrix. Among the 15-parameter family of realizations, a common choice is the canonic projection matrices P and P [START_REF] Luong | Canonic representations for the geometries of multiple projective views[END_REF]:

P ∼ (I (3×3) 0 (3×1) ) and P ∼ (H e ), (1)
where e is the second epipole defined as the right kernel of F, F T e ∼ 0 (3×1) and H ∼ [e ] × F is the canonic plane homography. This defines the canonic coordinate frame which is unique, provided normalization constraints for H and e . It will be seen in §2.3 that ||e || 2 = 1 and ||H || 2 = γ, where γ is an unknown constant scalar, is a convenient choice for my method. Note that ∼ means "equal up to scale" and [.] × is the cross-product skew-symmetric (3 × 3)-matrix. All entities are represented in homogeneous coordinates, i.e. are defined up to scale.

Relation to Previous Work

Most of the previous work on minimally parameterizing projective two-view motion deals with directly parameterizing the epipolar geometry. The fundamental matrix F is decomposed into the epipoles e and e and the epipolar transformation, which is a 1D projective transformation relating the epipolar pencils, represented by an homogeneous (2 × 2)-matrix g. Representing these entities with minimal parameters requires eliminating their homogeneity. This is usually done by normalizing each of them so that their largest entry is unity, which yields 3 possibilities for each epipole and 4 for the epipolar transformation, so 3 • 3 • 4 = 36 possible parameterizations.

Obtaining the fundamental matrix -or any other 2D entity such as the extended epipolar transformation [START_REF] Bartoli | Projective structure and motion from two views of a piecewise planar scene[END_REF] or the canonic plane homography H -from e, e and g requires then the use of 9 distinct parameterizations to model all cases [START_REF] Zhang | Determining the epipolar geometry and its uncertainty: A review[END_REF]. These cases coincide with 9 of the 36 previous ones.

In [START_REF] Luong | The fundamental matrix: Theory, algorithms and stability analysis[END_REF], the author proposes to restrict the two-view configurations considered to the cases where both epipoles are finite and can therefore be expressed in affine coordinates. Due to the homogeneity of the epipolar transformation, 4 distinct parameterizations are still necessary for g. A unique parameterization can then be used to form the fundamental matrix.

The method has been extended in [START_REF] Zhang | Determining the epipolar geometry and its uncertainty: A review[END_REF] to the general case, i.e. when the epipoles can be either finite or infinite. In this case, it is shown that the full 36 distinct parameterizations are necessary. This leads to a cumbersome and error-prone implementation of the optimization process.

In [START_REF] Bartoli | Projective structure and motion from two views of a piecewise planar scene[END_REF][START_REF] Zhang | Estimating the fundamental matrix by transforming image points in projective space[END_REF], the method has been revised so as to reduce the number of parameterizations using image transformations. In [START_REF] Bartoli | Projective structure and motion from two views of a piecewise planar scene[END_REF], the image transformations used are metric and the number of distinct parameterizations restricted to 3 plus one bilinear constraint on the entries of g, while in [START_REF] Zhang | Estimating the fundamental matrix by transforming image points in projective space[END_REF], the transformations used are projective, which allows one to reduce the number of parameterizations to 1. The main drawback is that in the transformed image space, the original noise model on the image features is not preserved. A means to preserve it, up to first order approximation, has been proposed in [START_REF] Zhang | Estimating the fundamental matrix by transforming image points in projective space[END_REF] for the gradient-weighted criterion, which is not the one used for bundle adjustment.

The Orthonormal Representation

Derivation. I introduce what I call the orthonormal representation of projective twoview motion. I consider the fundamental matrix representation of the motion. Any (3 × 3) rank-2 matrix is a fundamental matrix, i.e. represents a motion. Conversely, a projective two-view motion is represented by a unique fundamental matrix (up to scale). Therefore, deriving a minimal representation of projective two-view motion from its fundamental matrix representation implies considering two constraints; the rank-2 constraint and a normalization constraint, which fixes the relative scale of the fundamental matrix. Previous work has shown that these constraints are quite tricky to enforce directly on the fundamental matrix [START_REF] Bartoli | Projective structure and motion from two views of a piecewise planar scene[END_REF][START_REF] Luong | The fundamental matrix: Theory, algorithms and stability analysis[END_REF][START_REF] Zhang | Determining the epipolar geometry and its uncertainty: A review[END_REF][START_REF] Zhang | Estimating the fundamental matrix by transforming image points in projective space[END_REF].

To overcome this problem, instead of considering directly the fundamental matrix, I rather analyze its singular value decomposition F ∼ UΣV T , where U and V are (3 × 3) orthonormal matrices and Σ a diagonal one containing the singular values of F. The orthonormal representation is then derived while enforcing the constraints:

-rank-2: since F is a rank-2 matrix, Σ ∼ diag(σ 1 , σ 2 , 0) where σ 1 and σ 2 are strictly positive scalars and σ 1 ≥ σ 2 > 0; -normalization: since F is an homogeneous entity, I can scale its singular value decomposition such that

F ∼ U•diag(1, σ, 0)•V T where σ = σ 2 /σ 1 and 0 < σ ≤ 1 (σ 1 = 0 since F is rank-2).
This decomposition shows that any projective two-view motion can be represented by two (3 × 3) orthonormal matrices and a scalar. This gives the orthonormal representation of projective two-view motion as:

(U, V, σ) ∈ F where F ≡ O(3) 2 × {σ 0 < σ ≤ 1}, (2)
where O(3) is the Lie group of (3 × 3) orthonormal matrices. This representation is minimal in that it has 3 + 3 + 1 = 7 degrees of freedom. It can easily be computed from the singular value decomposition of the fundamental matrix. Note that σ = 1 may correspond to the case of an essential matrix, i.e. when cameras are calibrated. Any element of F represents a unique two-view motion since it can be used to recompose a unique fundamental matrix (see next paragraph), i.e. a (3 × 3) matrix where both the rank-2 and a normalization constraints have been enforced. However, a fundamental matrix has more than one orthonormal representations. For instance, given an orthonormal representation (U, V, σ) ∈ F, one can freely switch the signs of u 3 or v 3 while leaving the represented motion invariant. However, the space of fundamental matrices and the orthonormal representation are both 7-dimensional, which allows for minimal estimation.

Recovering 2D entities. The fundamental matrix corresponding to an orthonormal representation (U, V, σ) ∈ F can be recovered by simply recomposing the singular value decomposition:

F ∼ u 1 v T 1 + σu 2 v T 2 , (3)
where u i and v i are the columns of U and V respectively. Among all potential applications of the orthonormal representation, I will use it for bundle adjustment. Therefore, I will need to extract projection matrices from the fundamental matrix. This can be achieved directly from the orthonormal representation by recovering the second epipole and the canonic plane homography of [START_REF] Luong | Canonic representations for the geometries of multiple projective views[END_REF], equation [START_REF] Bartoli | Three new algorithms for projective bundle adjustment with minimum parameters[END_REF].

The second epipole is the last column of U: e ∼ u 3 , while the canonic plane homography can be formulated as

H ∼ [e ] × F ∼ [u 3 ] × u 1 v T 1 + σu 2 v T 2 . Since U is an O(3) matrix, [u 3 ] × u 1 = ±u 2 and [u 3 ] × u 2 = ∓u 1 which yields the canonic plane homography as H ∼ u 2 v T 1 -σu 1 v T 2
and the second projection matrix as:

P ∼ u 2 v T 1 -σu 1 v T 2 | u 3 . (4)
Normalization constraints discussed in §2.1 are clearly satisfied. This guarantees that the same canonic basis will be used through the optimization process.

Estimation Using Minimal Parameters

In this section, I use the previously described orthonormal representation of projective two-view motion to locally update a "base" estimate using the minimum 7 parameters.

Update using minimal parameters. Before going further, let us examine the case of 3D rotations. There does not seem to exist a minimal, complete and non-singular parameterization of the 3-dimensional set of rotations in 3D space. For example, consider the (3 × 3) rotation matrix R ∈ SO(3). Minimal representations, such as the 3 Euler angles θ lead to singularities. However, one can find representations that are locally minimal and non-singular, e.g. in a neighbourhood of R = I (3×3) , i.e. θ = 0 (3×1) . Therefore, most of the estimation processes of 3D rotations do not minimally parameterize the current estimate, but rather locally update its overparameterization. A typical example is to use a (3 × 3) rotation matrix representation that is locally updated as

R ← R • R(θ)
where R(θ) is any minimal and locally non-singular parameterization of 3D rotations, e.g. Euler angles. This method is used in [START_REF] Hartley | Euclidean reconstruction from uncalibrated views[END_REF] for the non-linear recovery of metric structure and motion, where R is updated and θ reset to zero after each iteration of the optimizer. The same scenario arises for the projective two-view motion. There does not seem to exist a minimal, complete and non-singular parameterization of the corresponding 7-dimensional space. Consequently, I propose locally updating a given estimate using minimal parameters.

Let us consider the orthonormal representation (U, V, σ) ∈ F of a given projective two-view motion. Each matrix (U, V) ∈ O(3) 2 can be locally updated using 3 parameters by considering the method described above for 3D rotations. The scalar σ ∈ {σ 0 < σ ≤ 1} is completely included into the parameterization. A means to ensure 0 < σ ≤ 1 is described below.

Completness. A first remark that immediately follows about the above-proposed method is whether all two-view configurations are covered. This arises from the fact that U and V are O(3) matrices, which may have positive or negative determinants, and are updated using SO(3) matrices, R(x) and R(y) respectively, which have only positive determinants. Therefore, the signs of det(U) and det(V) are not allowed to change during the optimization process. I claim that this is not a problem and that any fundamental matrix F can be reached from any initial guess F 0 , even if sign(det(U 0 )) = sign(det(U)), sign(det(V 0 )) = sign(det(V)) or both. To prove this claim, I show that any fundamental matrix F represented by (U, V, σ) has alternative representations where the signs of either det(U) or det(V), or both, have been switched. This is due to the non-uniqueness of the singular value decomposition, see e.g. [START_REF] Press | Numerical Recipes in C -The Art of Scientific Computing[END_REF]. Consider the recomposition equation ( 3) and observe that u 3 and v 3 , the third columns of U and V respectively, do not affect the result. Therefore, they are only constrained by the orthonormality of U and V. Hence, their signs can be arbitrarily switched, which accordingly switches the sign of the determinant of the corresponding matrix. For example, u 3 ← -u 3 switches the sign of det(U) while leaving invariant the represented fundamental matrix. This concludes the proof.

Implementation details. Through the iterations, σ may go outside of its boundaries. This is not a problem since the corresponding motion is still valid.

There are several possibilities to ensure the boundaries of σ such as using linear constraints. I propose to enforce these boundaries at each iteration while leaving the current estimate invariant. However, I have found during my experiments of § §5 and 6 that in practice, this does not affect the behaviour of the underlying optimization process. A way to proceed is:

-if σ < 0 then σ ← -σ, u 2 ← -u 2 , u 3 ← -u 3 ; -if σ > 1 then σ ← 1 σ , swap(u 1 ,u 2 ), u 3 ← -u 3 , swap(v 1 ,v 2 ), v 3 ← -v 3 .
One can easily check that these changes on the orthonormal representation leave the represented motion invariant.

Bundle Adjustment

In this section, I show how the orthonormal representation can be used for bundle adjustment of point features seen in two views. This is summarized in table 1. Similar results could be derived for other criteria, such as the minimization of the distances between points and epipolar lines or the gradient-weighted criterion [START_REF] Luong | The fundamental matrix: Theory, algorithms and stability analysis[END_REF][START_REF] Zhang | Determining the epipolar geometry and its uncertainty: A review[END_REF].

Two-view projective bundle adjustment expressed within the framework of [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], p.574 (algorithm A4.1). The initial guess of the fundamental matrix is denoted by F0.

Add the following steps: (i') Initialize the orthonormal representation (U, V, σ) by a scaled singular value decomposition of F0:

F0 ∼ U • diag(1, σ, 0) • V T .
(ii') Turn the full (r × 12) camera Jacobian matrix A = Ā into the minimal (r × 7) Jacobian matrix of the orthonormal representation:

A ← A • A ortho ,
where A ortho is given by equations [START_REF] Hartley | [END_REF][START_REF] Hartley | Euclidean reconstruction from uncalibrated views[END_REF];

Change the parameter update step as: (viii) Update the orthonormal representation as:

U ← U • R(x) V ← V • R(y) σ ← σ + δσ,
where δa T = (x T y T δσ) are the 7 motion update parameters, update the structure parameters by adding the incremental vector δ b and compute the new error vector;

Add the last step:

(xi) Return the computed F using equation (3) as:

F ∼ u1v T 1 + σu2v T 2 .
Table 1. Implementing my minimal estimator within the bundle adjustment Levenberg-Marquardt-based framework given in [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], p.574 (algorithm A4.1). Note that r is the number of residuals and that the second projection matrix have to be extracted from the orthonormal representation using equation (4) (for e.g. computing the error vector).

Cost function. Bundle adjustment consists in solving the following optimization problem, see e.g. [START_REF] Mclauchlan | Gauge invariance in projective 3D reconstruction[END_REF][START_REF] Triggs | Bundle ajustment -a modern synthesis[END_REF][START_REF] Zhang | Determining the epipolar geometry and its uncertainty: A review[END_REF]:

min a,b r T r,
where:

a and b are respectively motion and structure parameters (or parameters used to update them); -r is the vector of residual errors; -r 2 i = d 2 (q i , PQ i ) + d 2 (q i , P Q i ) is the i-th point residual error (d is the 2D Euclidean distance) corresponding to its reprojection error; -q i and q i are corresponding image points for the first and second images; -Q i are 3D reconstructed points and depend upon b; -P and P are projection matrices corresponding to the current motion estimate represented by a. They must correspond to a realization of the fundamental matrix. I have shown in §2.3, equation ( 4), how the canonic realization can be directly obtained from the orthonormal representation.

Analytical differenciation. Newton-type optimization methods, such as the widely used Levenberg-Marquardt, necessitate computing the Jacobian matrix J = (A | B) of the residual vector r with respect to motion and structure parameters a and b. While this can be achieved numerically using e.g. finite differences [START_REF] Press | Numerical Recipes in C -The Art of Scientific Computing[END_REF], it may be better to use an analytical form for both computational efficiency and numerical accuracy. I focus on the computation of A = ∂r ∂a since B = ∂r ∂b only depends upon structure parameterization. I decompose it as A (r×7) = Ā(r×12) • A ortho (12×7) where:

r is the number of residuals; -only the 12 entries of P are considered since P is fixed in the canonic reconstruction basis (1); -Ā = ∂r ∂p (p = vect(P ) where vect(.) is the row-wise vectorization) depends on the chosen realization of the fundamental matrices, i.e. on the coordinate frame employed. I have chosen the canonic projection matrices (1). This Jacobian matrix is employed directly for the overparameterization of [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. Deriving its analytical form is straightforward; -A ortho = ∂p ∂a is related to the orthonormal motion representation.

I therefore concentrate on deriving a closed-form expression for A ortho . If the minimal method of e.g. [START_REF] Zhang | Determining the epipolar geometry and its uncertainty: A review[END_REF] were used, 36 different Jacobian matrices corresponding to each parameterization would have to be derived. One of the advantages of my update scheme shown in table 1 and based on the orthonormal representation is that there exists a simple closed-form expression for A ortho .

Let us consider the orthonormal representation (U, V, σ). In this case, the motion update parameters are minimal and defined by a = (x 1 x 2 x 3 y 1 y 2 y 3 σ) T , where x = (x 1 x 2 x 3 )

T and y = (y 1 y 2 y 3 ) T are used to update U and V respectively. Since U and V are updated with respect to the current estimate, A ortho is evaluated at (U, V, σ), i.e. at a = a 0 = (0 (6×1) σ) T . Let U = U • R(x) and V = V • R(y) be the updated U and V. Equation ( 4) is used to derive a closed-form expression of the second canonic projection matrix after updating, i.e. corresponding to the orthonormal representation ( U, V, σ). By expanding, differentiating and evaluating this expression at a 0 , I obtain:

A ortho = ∂p ∂a = ∂p ∂(x 1 . . . y 3 σ) = ∂p ∂x 1 . . . ∂p ∂y 3 ∂p ∂σ , (5)
where:

∂p =                    vect(u 3 v 1 T | -u 2 ) • ∂x 1 vect(-σu 1 v 3 T | 0) • ∂x 2 vect(σu 3 v 2 T | u 1 ) • ∂x 3 vect(-u 2 v 3 T | 0) • ∂y 1 vect(-u 1 v 1 T -σu 2 v 2 T | 0) • ∂y 2 vect(u 2 v 2 T + σu 1 v 1 T | 0) • ∂y 3 vect(-u 1 v 2 T | 0) • ∂σ. (6)

The Multiple-View Case

In this section, I extend my projective two-view motion modelisation to multiple views. I analyse how to model additional views. I propose the orthonormal representation of projective multiple-view motion. As in the two-view case, this serves to devise elements for optimizing projective multiple-view motion over minimal parameters.

Modeling Additional Views

Once two views have been modeled, the coordinate frame is fixed. Therefore, an additional view does not have any free gauge and its complete projection matrix has to be modeled. Let P be such a (3 × 4) projection matrix. Since it is homogeneous, it has 11 degrees of freedom. This can be seen in several other ways. For example, one can consider the metric decomposition P ∼ K(R t). However, this decomposition is not available here since I deal with uncalibrated cameras. One may also interpret P ∼ (H αe), where H is a (3×3)-matrix, as a plane homography H with respect to the reference view, i.e. thus for which the projection matrix is (I (3×3) 0 (3×1) ), e a 3-vector that represents an epipole with the same view and α a scalar that accounts for the relative scale between H and e. Therefore, P has 8 + 2 + 1 = 11 degrees of freedom. This interpretation is related to that of plane+parallax, see e.g. [START_REF] Irani | Parallax geometry of pairs of points for 3d scene analysis[END_REF].

Relation to Previous Work

A common strategy for optimizing an homogeneous entity such as P is to overparameterize it by using all its entries. A normalization constraint is then softly imposed by using an hallucinated measurement, e.g. on the norm of P as ||P|| 2 -1 = 0. The drawback of this method is that more parameters than necessary are estimated, which increases the computational cost of the estimation process and may cause numerical instabilities. One could also renormalize P after each iteration as P ← P/||P|| 2 . Alternatively, one could fix one entry of P to a given value, e.g. 1, but this representation would have singularities.

The main drawback of these techniques is that a unique minimal parameterization does not suffice to express all cases. This leads to the necessity for multiple expressions of e.g. the Jacobian matrix for Newton-type optimizers, which might complexify implementation issues.

The Orthonormal Representation

The orthonormal representation of P ∼ (H αe) can be derived as follows. Let s = αe. This inhomogeneous 3-vector is a scaled version of e which has 3 degrees of freedom since it also encapsulates the relative scale α between H and e. Therefore, s ∈ R 3 and it can be directly parameterized by its 3 elements.

Consider now the homogeneous (3×3)-matrix H. As in the case of the fundamental matrix, see §2.3, I examine its singular value decomposition H ∼ UΣV T where U and V are (3 × 3) orthonormal matrices and Σ a diagonal one containing the singular values of H. Since H may be singular, see e.g. in §2.1 the canonic plane homography of [START_REF] Luong | Canonic representations for the geometries of multiple projective views[END_REF], but must not be rank-1 or null, Σ ∼ diag(1, σ 1 , σ 2 ), where 0 < σ 1 ≤ 1 and 0 ≤ σ 2 ≤ σ 1 . Therefore, the orthonormal representation of H writes as:

(U, V, σ 1 , σ 2 ) ∈ H where H ≡ F × {σ 2 0 ≤ σ 2 ≤ σ 1 }, and F ≡ O(3) 2 ×{σ 1 0 < σ 1 ≤ 1}, see equation (2).
As a byproduct, one can observe that I have derived the orthonormal representation H of 2D homographies which can be used as a starting point to devise minimal 8-parameter estimators for these transformations. Finally, I obtain the orthonormal representation of P, denoted by P as:

(U, V, σ 1 , σ 2 , s) ∈ P where P ≡ H × R 3 .
It is minimal in that it has 3 + 3 + 1 + 1 + 3 = 11 degrees of freedom.

Optimization With Minimal Parameters.

Mimicking the method of §2.4 for the projective two-view motion case, I obtain a way to minimally estimate projective multiple-view motion. Given a set of camera matrices, I represent two of them using the orthonormal projective two-view motion of §2.3. This fixes the coordinate frame. Each other view is then modeled by the orthonormal representation (U, V, σ 1 , σ 2 , s) ∈ P described above. Each O(3) matrix U, V can be updated using minimal parameters as e.g. U ← U • R(x) whereas σ 1 , σ 2 and s are directly optimized. As in the two-view case, one can derive algorithms to ensure the boundaries on σ 1 and σ 2 . A closed-form solution for the Jacobian matrix of the residuals with respect to the motion parameters can be derived in a manner similar to the twoview case. Another advantage of this representation is that one can directly compute the inverse of H, the (3 × 3) leading part of P, from its orthonormal representation. This may be useful for e.g. projecting 3D lines or estimating 2D homographies using a symmetric cost function.

Experimental Results Using Simulated Data

In this section, I compare the algorithm using the orthonormal motion representation (see table 1 and § §2.4 and 4.4) to existing ones using simulated data. I use the Levenberg-Marquardt method to perform optimization. Points are minimally parameterized as in [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], p.579. The test bench consists of 50 points lying inside a sphere with a radius of 1 meter observed by cameras with a focal length of 1000 (expressed in number of pixels). Each of these cameras is looking at the center of the sphere and is situated at a distance of 10 meters from it. The baseline between two consecutive cameras is 1 meter.

Points are generated in 3D space, projected onto the images and corrupted by an additive centered gaussian noise with varying standard deviation.

I measure the two quantities characteristic of a bundle adjustment process, computational cost, i.e. CPU time to convergence and the error at convergence, versus the standard deviation of added image noise. I also measure the error of the current estimate as a function of time through the optimization processes. The plots correspond to median values over 300 trials. The bundle adjustments are initialized by the values obtained using the 8 point algorithm [START_REF] Hartley | In defence of the 8-point algorithm[END_REF] and the triangulation method described in [START_REF] Hartley | [END_REF] for two views. Each other view is then registered in turn by linearly computing its camera matrix. I compare the following algorithms:

-ORTHO: uses the minimal methods given in table 1 and § §2.4 and 4.4; -MAPS and IMROT (Image Rotation): are other minimal methods given in [START_REF] Bartoli | Three new algorithms for projective bundle adjustment with minimum parameters[END_REF], the associated research report. These methods are equivalent to those described in [START_REF] Zhang | Determining the epipolar geometry and its uncertainty: A review[END_REF][START_REF] Zhang | Estimating the fundamental matrix by transforming image points in projective space[END_REF] in the sense that the number of unknowns is minimal; -FREE: uses an overparameterization with free gauges, namely all the entries of the camera matrices are optimized; -NORMALIZED: uses an overparameterization plus hallucinated measurements to prevent the gauge to drift [START_REF] Mclauchlan | Gauge invariance in projective 3D reconstruction[END_REF]; -PARFREE (Partially Free): uses a partially free gauge by completely parameterizing all camera matrices except the first one [START_REF] Hartley | Projective reconstruction and invariants from multiple images[END_REF]; I conduct a serie of experiments using two cameras. One can observe on figure 1 (b) that, roughly speaking, all methods converge to the same minimum of the cost function.

Methods that have a slightly less reliable convergence than the others are FREE and PARFREE.

Figure 1 (a) shows that, for roughly the same convergence properties, there are quite big discrepancies between the computational cost of each method. The method that has the highest computational cost is FREE, followed by NORMALIZED and PARFREE. This can be explained by the fact that these methods have more unknowns to estimate than the minimal ones. This requires more computational time for each iteration to be performed. Finally, methods using the minimal number of parameters, MAPS, ORTHO and IMROT have the lowest computational cost, roughly the same.

In the light of these results, it is clear that methods using minimal parameters should be preferred for both computational cost and convergence properties. The method OR-THO, relying on the orthonormal representation given in this paper has the advantage of simplicity. However, in order to understand and explain the behaviour of the different methods, I have measured the number of iterations and the computational cost of these iterations. These results are shown on figure 2. In more detail, I have found that methods FREE or PARFREE, leaving the gauge drift freely have very bad convergence properties, performing more iterations, roughly twice more, than the others, see figure 2 (a). Method NORMALIZED performs a num-ber of iterations smaller than all the other methods but involves solving a much more costly linear system at each iteration, see table 2. Methods using the minimal number of parameters are trade-offs between the number of iterations and their computational cost: each iteration has a low computational cost and the number of iterations needed is in-between those of free gauge methods and NORMALIZED. This explains why these methods achieve the lowest total computational cost.

Figure 2 (b) shows the evolution of reprojection error for the different optimization processes. This experiment is useful in the sense that the time to convergence previously measured is highly dependent on how convergence is determined, e.g. by thresholding two consecutive errors, and does not account for the ability of the algorithms to quickly, or not, reach almost the value of convergence. This experiment has been conducted using the same test bench as previously with a noise level on image point positions of 0.5 pixel. One can see on this figure that methods based on a minimal parameterization reach their value of convergence before the others. The NORMALIZED and PARFREE methods take roughly twice the same time, while the FREE method takes three times more.

Finally, I conduct experiments using 10 views. I observe that the differences between the algorithms observed in the two-view case are decreased while those requiring the lowest computation time are the same, i.e. MAPS, ORTHO and IMROT. Other experiments on the convergence properties of the algorithms in the multiple-view case yield conclusions similar to the two-view case.

Experiments on Real Images

In this section, I validate my algorithms using real images. I first consider the case of two images. In order to cover all possibilities for the epipoles to be close to the images or at infinity, I use pairs of the images shown on figure 3. Initial values for structure and motion are computed as in the case of simulated data. Results are shown in table 3. For each combination of image pair and each algorithm, I estimate the CPU time to convergence T and the error at convergence E. The last row of the table show the mean values T and Ē of T and E for each algorithm over the set of image pairs. These results confirmed those obtained using simulated data. I have also tested the algorithms on all four images of figure 3. Initial values have been obtained by registering each view to an initial guess of structure and motion obtained from the two first ones. The results are the followings: all algorithms converge with a final error of 0.73 pixels and their relative performances in terms of computation times to convergence were equivalent to those obtained in the case of two views.

Conclusions

I studied the problem of optimizing projective motion over minimal sets of parameters. I proposed the orthonormal representation of projective two-view motion. I showed how this can be used to locally update projective two-view motion using a minimal set of 7 parameters. The canonic projection matrices can be directly extracted from the orthonormal representation. I extent this representation to projective multiple-view motion. As a byproduct, I derive the orthonormal representation of 2D homographies. The method can be plugged into most of the (possibly sparse) non-linear optimizers such as Levenberg-Marquardt. I gave a closed-form expression for the Jacobian matrix of the residuals with respect to the motion parameters, necessary for Newton-type optimization techniques.

The introduced orthonormal representation seems to be a powerful tool for minimal optimization of homogeneous entities in particular.

I conducted experiments on simulated and real data. My conclusions are that methods based on minimal parameter sets perform better than the others, in terms of computational cost while achieving equivalent results in terms of convergence properties. The most interesting results are obtained in the two-view case. Existing algorithms that do not constrain the gauge by any means perform worse than the others.

I will make an implementation of the method available on my home-page. In future work, I plan to investigate the use of the orthonormal representation introduced in this paper to model other algebraic entities and devise minimal estimation techniques for them.

Fig. 1 .

 1 Fig. 1. Comparison of (a): the CPU time to convergence and (b): value of the error at convergence versus varying image noise for different methods.

Fig. 2 .

 2 Fig. 2. Comparison of (a): the number of iterations to convergence versus varying image noise and (b): the evolution of reprojection errors.

Fig. 3 .

 3 Fig. 3. Real images used to validate the algorithms.

Table 2 .

 2 CPU time per iteration (second) for each method.

	MAPS ORTHO IMROT FREE NORMALIZED PARFREE
	0.3841 0.3838 0.3864 0.4741	0.7043	0.4717

Table 3 .

 3 A, B 0.63 2.45 0.63 2.39 0.63 2.47 0.68 3.98 0.63 2.99 0.68 3.02 A, C 0.71 2.38 0.71 2.41 0.71 2.40 0.77 4.01 0.71 3.56 0.71 3.71 -∞ -∞ A, D 0.45 2.03 0.45 1.76 0.45 2.19 0.57 3.13 0.45 3.09 0.45 2.93 ∞ ∞ B, C 0.88 3.53 0.88 3.39 0.88 3.55 1.23 6.70 0.88 5.12 0.88 4.63 ∞ -∞ B, D 0.59 2.33 0.59 2.10 0.59 2.81 0.59 3.99 0.59 3.41 0.59 3.56 C, B 0.51 1.91 0.51 1.92 0.51 2.02 0.51 3.39 0.51 2.79 0.51 3.04 average Ē and T 0.628 2.430 0.628 2.328 0.628 2.573 0.725 4.200 0.628 3.493 0.637 3.482 Error at convergence E and time to convergence T obtained when combining pairs of images from figure 3.

	epipoles views e e	MAPS E T	ORTHO E T	IMROT E T	FREE E T	NORMALIZED PARFREE E T E T
	-∞ ∞


Failure of the raw ensemble Kalman filter (EnKF) ◮ With the exception of Gaussian and linear systems, EnKF fails to provide a proper estimation on most systems.

◮ To properly work, it needs fixes: localisation and inflation. ◮ EnKF relies for its analysis on the first and second-order empirical moments:

x = 1 N N ∑ n=1 x k , P = 1 N -1 N ∑ k=1 (x k -x) (x k -x) T .
(1)

Yet, x and P may not be the true moments of the true filtering distribution (assuming there is one!). Hidden true moments of the true filtering distribution: x b and B.

M. Bocquet Int. Conf. Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012 3 / 29

Getting more from the ensemble ◮ Idea: even under Gaussian assumptions of the true distribution, the pdf p(x|x 1 , . . . , x N ) extracts more information than p(x|x, P).

◮ Using Gaussian assumptions, and being only interested in the filtering problem, one can get (hierarchical reasoning):

p(x|x 1 , . . . , x N ) = 1 p(x 1 , . . . , x N ) dx b dB p(x|x b , B)p(x 1 , . . . , x N |x b , B)p(x b , B) .

(2)

◮ p(x|x b , B): the standard Gaussian prior but based on the true statistics. Choosing priors for the background statistics ◮ To progress, we need to make assumptions on the background statistics p(x b , B): the statistics of the error statistics or hyperpriors.

◮ p(x 1 , . . . , x N |x b , B) = ∏ N k=1 p(x k |x b , B) ◮ p(x b ,
A very simple choice is a weakly informative prior: the Jeffreys' prior [Jeffreys 1961] with an additional assumption of independence for x b and B:

p(x b , B) ≡ p J (x b , B) = p J (x b )p J (B) and p J (x b ) = 1 , p J (B) = |B| -M+1 2 .
◮ With Jeffreys prior, it is possible to perform the integral (with additional complications due to rank-deficiency usually not dealt with by mathematicians). Principle of the EnKF-N ◮ The prior of EnKF and the prior of EnKF-N:

p(x|x, P) ∝ exp - 1 2 (x -x) T P -1 (x -x) p(x|x 1 , x 2 , . . . , x N ) ∝ (x -x) (x -x) T + ε N (N -1)P -N 2 , (3)
with ε N = 1 (mean-trusting variant), or ε N = 1 + 1 N (original variant).

◮ Ensemble space decomposition (ETKF version of the filters): x = x + Aw.

◮ The variational principle of the analysis (in ensemble space): Compute Y = HA, δ = y -Hx

J (w) = 1 2 (y -H(x + Aw)) T R -1 (y -H(x + Aw)) + N -1 2 w T w J (w) = 1 2 (y -H(x + Aw)) T R -1 (y -H(x + Aw)) + N 2 ln ε N + w T w . (4)
Find the minimum:

w a = min w (δ -Yw) T R -1 (δ -Yw) + N ln ε N + w T w Compute x a = x + Aw a .
Compute The primal EnKF-N Forced 2D turbulence model

Ω a = Y T R -1 Y + N (εN+w T a w a) I N -2w a w T a (ε N +w T a w a ) 2 -1 Compute W a = {(N -1)Ω a } 1/2 U Compute x a k = x a + A W a k M.
dx m dt = (x m+1 -x m-2 )x m-1 -x m + F ,
◮ Forced 2D turbulence model ∂ q ∂ t + J(q, ψ) = λ q + ν∆ 2 q + F , q = ∆ψ , (5)
where J(q, ψ) = ∂ x q∂ y ψ -∂ y q∂ x ψ, q is the vorticity 2D field, ψ is the current function 2D field, F is the forcing, λ amplitude of the friction, ν amplitude of the biharmonic diffusion, grid: 64 × 64 small enough to be in the sufficient-rank regime.

◮ Setup of the experiment: Time-lag between update: ∆ t = 2, fully observed, R = 0.1I. Kuramato-Sivashinski model, [Kuramato, 1975;Sivashinski, 1977] NEDyM economical model, Lagrangian duality ◮ The primal EnKF-N cost function:

J (w) = 1 2 (y -H(x + Xw)) T R -1 (y -H(x + Xw)) + N 2 ln ε N + w T w . (6)
◮ Idea: Split the radial degree of freedom of w, that is √ w T w, from its angular degrees of freedom, that is w/ √ w T w.

◮ Lagrangian:

L (w, ρ, ζ ) = 1 2 (δ -Yw) T R -1 (δ -Yw) + 1 2 ζ w T w -ρ + N 2 ln (ε N + ρ) , (7)
where δ = y -Hx.

◮ Saddle point equations:

ζ ⋆ = N/(ε N + ρ ⋆ ) ζ ⋆ w ⋆ = -Y T R -1 (δ -Yw ⋆ ) ⇒ ρ ⋆ = N/ζ ⋆ -ε N w ⋆ = ζ ⋆ + Y T R -1 Y -1 Y T R -1 δ (8) M.
D(ζ ) = inf w sup ρ≥0 L (w, ρ, ζ ) = 1 2 δ T R + Yζ -1 Y T -1 δ + ε N ζ 2 + N 2 ln N ζ - N 2 . ( 9
)
◮ Dual and primal problems:

∆ = inf ζ >0 D(ζ ) and Π = inf w J (w) . (10)
◮ Strong duality result (non quadratic, non-convex case!): The dual EnKF-N scheme

∆ = Π . (11)
Requires: The forecast ensemble {x k } k=1,...,N , the observations y, and error covariance matrix R

Compute the mean x and the anomalies A from {x k } k=1,...,N .

Compute Y = HA, δ = y -Hx

Find the minimum:

ζ a = min ζ ∈]0,N/ε N ] δ T R + Yζ -1 Y T -1 δ + ε N ζ + N ln N ζ -N (12)
Compute

w a = Y T R -1 Y + ζ a -1 Y T R -1 δ . Compute x a = x + Aw a .
Compute ◮ The iterative extended Kalman filter [Wishner et al., 1969;Jazwinski, 1970] IEKF ◮ The iterative extended Kalman smoother [Bell, 1994] IEKS ✞ ✝

Ω a = Y T R -1 Y + ζ a 2ε N N ζ a -1 -1 Compute W a = {(N -1)Ω a } 1/2 U Compute x a k = x a + AW a k M.

☎ ✆

Much too costly + needs the TLM and the adjoint -→ ensemble methods ◮ The iterative ensemble Kalman filter [START_REF] Sakov | An iterative EnKF for strongly nonlinear systems[END_REF][START_REF] Bocquet | Combining inflation-free and iterative ensemble kalman filters for strongly nonlinear systems[END_REF] ◮ The iterative extended Kalman filter [Wishner et al., 1969;Jazwinski, 1970] IEKF ◮ The iterative extended Kalman smoother [Bell, 1994] IEKS ✞ ✝

☎ ✆

Much too costly + needs the TLM and the adjoint -→ ensemble methods ◮ The iterative ensemble Kalman filter [START_REF] Sakov | An iterative EnKF for strongly nonlinear systems[END_REF][START_REF] Bocquet | Combining inflation-free and iterative ensemble kalman filters for strongly nonlinear systems[END_REF] ◮ The iterative extended Kalman filter [Wishner et al., 1969;Jazwinski, 1970] IEKF ◮ The iterative extended Kalman smoother [Bell, 1994] IEKS ✞ ✝

☎ ✆

Much too costly + needs the TLM and the adjoint -→ ensemble methods ◮ The iterative ensemble Kalman filter [START_REF] Sakov | An iterative EnKF for strongly nonlinear systems[END_REF][START_REF] Bocquet | Combining inflation-free and iterative ensemble kalman filters for strongly nonlinear systems[END_REF] Iterative ensemble Kalman filters ◮ The iterative extended Kalman filter [Wishner et al., 1969;Jazwinski, 1970] IEKF ◮ The iterative extended Kalman smoother [Bell, 1994] IEKS ✞ ✝

☎ ✆

Much too costly + needs the TLM and the adjoint -→ ensemble methods ◮ The iterative ensemble Kalman filter [START_REF] Sakov | An iterative EnKF for strongly nonlinear systems[END_REF][START_REF] Bocquet | Combining inflation-free and iterative ensemble kalman filters for strongly nonlinear systems[END_REF] Iterative ensemble Kalman filters ◮ A fairly recent idea: [START_REF] Gu | An iterative ensemble Kalman filter for multiphase fluid flow data assimilation[END_REF]: The idea.

[Kalnay & Yang, 2010]: A step in the right direction. [Sakov, Oliver & Bertino, 2011]: The "pièce de résistance" [START_REF] Bocquet | Combining inflation-free and iterative ensemble kalman filters for strongly nonlinear systems[END_REF]: Correction of the bundle scheme + ensemble transform form.

◮ IEnKF cost function in ensemble space:

J (w) = 1 2 (y 2 -H 2 (M 1→2 (x 1 + A 1 w))) T R -1 2 (y 2 -H 2 (M 1→2 (x 1 + A 1 w))) + 1 2 (N -1)w T w . (13)
◮ Gauss-Newton scheme:

w (p+1) = w (p) -H -1 (p) ∇ J (w (p) ) ∇ J (p) = -Y T (p) R -1 2 y 2 -H 2 M 1→2 (x + A 1 w (p) ) + (N -1)w (p) , H (p) = (N -1)I N + Y T (p) R -1 2 Y (p) , Y (p) = [H 2 M 2←1 A 1 ] ′ (p) , (14)

◮

  where F = 8, and the boundary is cyclic. Chaotic dynamics, topological dimension of 13, a doubling time of about 0.42 time units, and a Kaplan-Yorke dimension of about 27.1. ◮ Setup of the experiment: Time-lag between update: ∆ t = 0.05 (about 6 hours for a global model), fully observed, R = I. . Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012 8 / 29 Application to the Lorenz '95 model ◮ EnKF-N: analysis rmse versus ensemble size, for ∆t = 0.05. Local version: LETKF-N, with N = 10 (beware ∆t = 0.01 requires a correction). (optimal localisation) LETKF (optimal inflation and localization) M. Bocquet Int. Conf. Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012 10 / 29

◮◮

  Comparison of: EnKF with uniform inflation, EnKF-N, adaptive inflation EnKF (EnKF-ML), N = 80 (rank-sufficient regime). Starting away or close from the truth. EnKF-N also tested on:Lorenz '63 model,[Lorenz, 1963]

  ://images.math.cnrs.fr M. Bocquet Int. Conf. Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012 13 / 29

Finite-size iterative

  ensemble Kalman smoothers ◮ Setup: Lorenz '95, M = 40, N = 20, ∆t = 0.30, R = I. ◮ Comparison of EnKF-N, IEnKF-N, IEnKS-N, ETKS-N, with L = 10. ◮ Lin-IEnKS-N has (understandably) diverged. . Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012 25 / 29 The iterative ensemble Kalman filter & smoother Finite-size iterative ensemble Kalman smoothers ◮ Setup: 2D turbulence, 64 × 64, N = 40, ∆t = 2, R = 0.1I. ◮ Comparison of EnKF-N, IEnKF-N, IEnKS-N, ETKS-N, with L = 20.

  Compute the mean x and the anomalies A from {x k } k=1,...,N .

	EnKF-N: algorithm
	Requires: The forecast ensemble {x k } k=1,...,N , the observations y, and error covariance matrix R
	M. Bocquet	Int. Conf. Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012	6 / 29

  The iterative extended Kalman filter[Wishner et al., 1969; Jazwinski, 1970] IEKF ◮ The iterative extended Kalman smoother[Bell, 1994] IEKS

	✞	☎
	✝	✆

Bocquet Int. Conf. Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012 Iterative ensemble Kalman filters ◮ Much too costly + needs the TLM and the adjoint -→ ensemble methods M. Bocquet Int. Conf. Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012

  IEnKF ◮ The iterative ensemble Kalman smoother [This talk. . . ] IEnKS

	✞ ✝ It's TLM and adjoint free!	☎ ✆

M. Bocquet Int. Conf. Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012 19 / 29

  IEnKF ◮ The iterative ensemble Kalman smoother [This talk. . . ] IEnKS

	✞ ✝ It's TLM and adjoint free!	☎ ✆
	✞ ✝ Don't want to be bothered by inflation tuning?	☎ ✆

M. Bocquet Int. Conf. Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012 19 / 29

  IEnKF ◮ The iterative ensemble Kalman smoother [This talk. . . ] IEnKS The finite-size iterative ensemble Kalman filter [Bocquet and Sakov, 2012] IEnKF-N ◮ The finite-size iterative ensemble Kalman smoother [This talk. . . ] IEnKS-N

	✞ ✝ It's TLM and adjoint free!	☎ ✆
	✞ ✝ Don't want to be bothered by inflation tuning?	☎ ✆
	◮

M. Bocquet Int. Conf. Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012 19 / 29

  Finite-size iterative ensemble Kalman filters ◮ Setup: Lorenz '95, M = 40, N = 40, ∆t = 0.05 -0.60, R = I. Finite-size iterative ensemble Kalman smoothers ◮ Setup: Lorenz '95, M = 40, N = 20, ∆t = 0.05, R = I.

	◮ Comparison of EnKF-N, EnKF (optimal inflation), IEnKF-N (bundle and transform), ◮ Comparison of EnKF-N, IEnKS-N, Lin-IEnKS-N, EnKS-N, IEnKS-N with a large
	IEnKF (bundle and transform, optimal inflation) inflation (reduced-rank 4D-Var?), with L = 20.
		0.7 0.9 0.22 0.8 0.20			N=40			EnKF-N Lin-IEnKS-N EnKS IEnKS-N λ=10
								IEnKS-N
		0.6 0.18
	Re-analysis rmse	0.3 0.4 0.5 0.16 Analysis rmse 0.14 0.12						EnKF-N
								EnKF (opt. infl.)
								IEnKF (bundle, opt. infl.)
		0.10						IEnKF-N (transform)
		0.2 0.08						IEnKF (transform, opt. infl.) IEnKF-N (bundle)
		0.15	0	0.1		0.2	0.3	0.4	0.5	0.6
			0	2	4	6 Time interval between updates 8 10 12 14 Lag (number of cycles)	16	18	20

M. Bocquet Int. Conf. Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012 20 / 29 M. Bocquet Int. Conf. Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012 22 / 29 M. Bocquet Int. Conf. Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012 24 / 29

Outline

The primal EnKF-N Iterative ensemble Kalman filters ◮ The iterative extended Kalman filter [Wishner et al., 1969;Jazwinski, 1970] ◮ Finite-size versions of the filter are just defined by substituting the prior:

◮ As a variational reduced method, one can use Gauss-Newton [START_REF] Sakov | An iterative EnKF for strongly nonlinear systems[END_REF],

Levenberg-Marquardt [START_REF] Bocquet | Combining inflation-free and iterative ensemble kalman filters for strongly nonlinear systems[END_REF][START_REF] Chen | Ensemble randomized maximum likelihood method as an iterative ensemble smoother[END_REF], etc, minimisation schemes (not limited to quasi-Newton).

◮ Essentially a lag-one smoother. Does the job of a lag-one 4D-Var, with dynamical error covariance matrix and without the use of the TLM and adjoint! Very efficient in very nonlinear conditions if one can afford the multiple ensemble propagations.

M. Bocquet Int. Conf. Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012

Iterative ensemble Kalman smoothers ◮ In a mildly nonlinear context (built on linear and Gaussian hypotheses) Many earlier studies, see [START_REF] Cosme | Smoothing problems in a bayesian framework and their linear gaussian solutions[END_REF] for a review, and [START_REF] Cosme | Implementation of a reduced-rank, square-root smoother for ocean data assimilation[END_REF] for an application to oceanography.

◮ In a non-sequential but very non-linear context

Many earlier studies, for instance [START_REF] Evensen | An ensemble Kalman smoother for nonlinear dynamics[END_REF]van Leeuwen, 2000] [Chen &Oliver, 2012] in the context of reservoir modelling

The IEnKS cost function is just the extension of the IEnKF cost function for a temporal window of L cycles.

M. Bocquet

Int. Conf. Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012

Conclusions

A new prior for the ensemble forecast meant to be used in an EnKF analysis has been built. It takes into account sampling errors.

It yields a new class of filters EnKF-N, that does not seem to require inflation supposed to account for sampling errors.

Local variants (both LA and CL) available.

Dual variant EnKF-N is an EnKF with built-in optimal inflation (accounting for sampling errors).

Almost linear regime more problematic because of Jeffreys' prior. Another hyperprior is needed.

The iterative ensemble Kalman filter has been generalised to an iterative ensemble Kalman smoother (IEnKF). It is an En-Var method.

It is tangent linear and adjoint free. It is, by construction, flow-dependent.

Though based on Gaussian assumptions, it can offer better retrospective analysis than standard Kalman smoothers in midly nonlinear conditions.

When affordable, it beats other Kalman filter/smoothers in strongly non-linear conditions.

M. Bocquet Int. Conf. Ens. Meth. Geo. Sc., Toulouse, France, 12-16 November 2012 28 / 29


Marc-André Ameller

Adrien Bartoli

Long Quan

Minimal Metric Structure and Motion From Three Affine Images

Keywords:

Structure and motion from minimal data is essential to bootstrap robust methods based on random sampling such as RANSAC or LMS. Let us consider the affine camera model and make the hypotheses of zero skew and unity aspect ratio. In this case, at least 4 points in 3 images are necessary to recover structure and motion. We propose a parametrization based on metric structure rather than camera motion parameters which have been previously used. The structure of 4 points is computed in closed-form by solving a quadratic equation. Unstable configurations are also investigated. Experimental results on simulated data and real images demonstrate that our algorithm allows fast estimation when included in a robust estimation process.

Introduction

Obtaining 3D structure and motion from image feature correspondences only is a fundamental task in computer vision. Particularly, reconstructing 3D points or estimating relative camera motions from minimal data is of primary importance for various numerical estimation procedures such as robust algorithms. In this paper, we propose an algorithm to obtain metric reconstruction from three affine images with the minimum of four point correspondences. While structure and motion from two affine views is ambiguous, it becomes possible when three or more are available. For that reason, this topic has been studied in the case of point [START_REF] Ullman | The Interpretation of Visual Motion[END_REF][START_REF] Tomasi | Shape and motion from image streams under orthography: A factorization method[END_REF] or line [START_REF] Astrom | Structure and motion from lines under affine projections[END_REF] features. A calibrated affine camera may be modeled by either orthographic, weak perspective or para-perspective projection [START_REF] Poelman | A paraperspective factorization method for shape and motion recovery[END_REF].

In this paper, a polynomial formulation of this minimal reconstruction problem is proposed based on a direct parametrization of the metric structure of a set of four points. Using tools from algebraic geometry, a closed-form solution to reconstruction and unstable cases are all given in the same framework. The main advantage of the method is that it gives directly stable structure by solving a quadratic equation rather than the less well-defined motion parameters.

The paper is organized as follow. Section 2 summarizes the affine camera model. Section 3 formulates the problem in terms of a polynomial system. Experimental results on simulated and real data are respectively given in sections 4 and 5.

The affine camera model

The general perspective camera model is described by a 3 × 4 matrix. Consider a point M ∼ (X, Y, Z, T ) T of the projective space, expressed in homogeneous coordinates. Its projection m ∼ (x, y, z) T on the image (considered as a projective plane) is given by the formula:

  x y z   ∼   p 1,1 p 1,2 p 1,3 p 1,4 p 2,1 p 2,2 p 2,3 p 2,4 p 3,1 p 3,2 p 3,3 p 3,4       X Y Z T     ,
where ∼ means "equal up to scale". In this paper, we deal with the affine camera model. This means that points at infinity in 3D space are projected to points at infinity in the image. This implies that:

p 3,1 = p 3,2 = p 3,3 = 0.
The projection equation can then be reduced to:

1 z x y = p 1,1 p 1,2 p 1,3 p 2,1 p 2,2 p 2,3   X/T Y /T Z/T   + t 1 t 2 .
Note that this is now a strict equality. Using decompostion QR [START_REF] Quan | Self-calibration of an affine camera from multiple views[END_REF]:

u v = λ 1 0 ξ 1 0 0 0 1 0 R   X Y Z   + t 1 t 2 , ( 1
)
where R is a rotation, λ is a scale factor, ξ is the aspect ratio, and is the skew factor. Camera calibration is equivalent to the knowledge of λ, ξ and . If and ξ are known, we can perform a change of image coordinates so that the projection writes:

u v = λ R T 1 R T 2   X Y Z   , (2)
where R T 1 and R T 2 are two rows of a 3 × 3 rotation matrix. In practice, we assume = 0 and ξ = 1. These assumptions are valid for most modern cameras.

Metric reconstruction

Problem formulation

Let us consider three affine cameras, each one represented by its projection matrix:

P i = λ i R i 1 T R i 2 T
, where i ∈ {1, 2, 3}. For each camera i, we can define an arbritrary plane P i , which represents its focal plane (the plane where 3D points are orthogonally projected). These planes are defined up to a translation. We expressed the recovered structure and motion in the first camera frame, which implies

P 1 = (I 2 |0 2 ).
The minimal configuration. Let us count the number of degrees of freedom of the problem. The first camera may be fixed arbitrarily, and for each additionnal camera, 6 additionnal d.o.f. corresponding to a 3 d.o.f. rotation and a 2 d.o.f. translation (translation is meaningful only parallel to the image plane), and a 1 d.o.f. scale factor. Next, we consider the number of constraints that might be given by each corresponding points: each point gives 6 constraints and 3 unknowns. The problem with 3 cameras and n points may have a solution if and only if

3n + 6 × (3 -1) ≤ 6 × n.
This implies that we need at least 4 points to solve the problem with 3 cameras. This is the minimal configuration.

Solving the Problem

The fundamental equation. Now, for each point M in space, we denote δ i (M ), the distance between M and P i . Then for each pair (M p , M q ) of points in space, and for the i-th camera, using the Pythagore theorem, we have:

λ i ||m -m || q p i i λ i ||m -m || q p i i m i p i q m M p δ (Μ ) i p δ (Μ ) q i ||Μ Μ || p q P i M q Figure 1. Deriving equation (3) |δ i (M p ) -δ i (M q )| 2 + λ 2 i m i p -m i q 2 = M p M q 2 ,( 3
)
where m i p = P i M p , and m i q = P i M q (see figure 1). This equation have interesting properties.

For our reconstruction problem, the calibration of the camera is known so that P i could be written as in [START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF]. The coordinates of corresponding points in images are known. Only the distances m i p -m i q are known. Eliminating scale factors. It is clear that the problem is defined up to a global scaling factor, i.e. multiplying all distances by any non-zero constant does not change the 3D structure, λ 1 can therefore be fixed to unity. Using 4 corresponding points, the affine fundamental matrix can be computed [START_REF] Shapiro | 3D motion recovery via affine epipolar geometry[END_REF] between cameras 1 and i, for each i = 1. The scale factors λ i can also be computed from the affine fundamental matrices as given in [START_REF] Shapiro | 3D motion recovery via affine epipolar geometry[END_REF].

Reduction of the fundamental equation. For simplicity, we introduce the following notations :

α i,p,q = λ 2 i m i p -m i q 2 X p = δ 1 (M p ) Y p = δ 2 (M p ) Z p = δ 3 (M p ).
Equation ( 3) is then rewritten as:

(X p -X q ) 2 + α 1,p,q = M p M q 2 (Y p -Y q ) 2 + α 2,p,q = M p M q 2 (Z p -Z q ) 2 + α 3,p,q = M p M q 2 .
Writing the above equations for cameras i and j, (i = j), and substracting them, we obtain:

(X p -X q ) 2 + α 1,p,q = (Y p -Y q ) 2 + α 2,p,q (X p -X q ) 2 + α 1,p,q = (Z p -Z q ) 2 + α 3,p,q .
Remembering that X p , Y p , Z p are defined up to the choice of an arbitrary plane P i , then up to a translation, we can set X 0 = Y 0 = Z 0 = 0. Consider a system given by 4 points in 3 images, it is described by the following polynomials:

P 1,2 = X 2 1 -Y 2 1 + k 1,2 P 1,3 = X 2 2 -Y 2 2 + k 1,3 P 1,4 = X 2 3 -Y 2 3 + k 1,4 P 2,3 = (X 1 -X 2 ) 2 -(Y 1 -Y 2 ) 2 + k 2,3 P 2,4 = (X 1 -X 3 ) 2 -(Y 1 -Y 3 ) 2 + k 2,4 P 3,4 = (X 2 -X 3 ) 2 -(Y 2 -Y 3 ) 2 + k 3,4 Q 1,2 = X 2 1 -Z 2 1 + l 1,2 Q 1,3 = X 2 2 -Z 2 2 + l 1,3 Q 1,4 = X 2 3 -Z 2 3 + l 1,4 Q 2,3 = (X 1 -X 2 ) 2 -(Z 1 -Z 2 ) 2 + l 2,3 Q 2,4 = (X 1 -X 3 ) 2 -(Z 1 -Z 3 ) 2 + l 2,4 Q 3,4 = (X 2 -X 3 ) 2 -(Z 2 -Z 3 ) 2 + l 3,4 , (4)
where:

k i,j = α 1,i-1,j-1 -α 2,i-1,j-1 , l i,j = α 1,i-1,j-1 -α 3,i-1,j-1 .
The goal is to compute X i , Y j and Z k from k i,j and l i,j . We first look at unstable cases for which small errors on the input data may give unreasonably large errors on the solution. We finally give a complete algorithm for structure and motion recovery.

Instability Conditions

To study the unstable cases, we introduce a new polynomial system obtained with pertubated measurement coefficients k i,j and l i,j :

k meas i,j = k i,j + δk i,j l meas i,j = l i,j + δl i,j .
Let the new solution corresponding to these pertubated measurements be:

X meas i = X i + δX i Y meas i = Y i + δY i Z meas i = Z i + δZ i .
Assume X i , Y i , Z i are the true solution of the original system, after a few algebraic manipulations, the pertubated system can be reduced in matrix form to:

δK = NδX + ε where δK =         δk 1,2 δk 1,3 δk 2,3 δl 1,2 δl 1,3 δl 2,3         and δX =         δX 1 δX 2 δY 1 δY 2 δZ 1 δZ 2         .
The matrix N has entries linear in X i , Y j , Z k and ε is a vector whose coefficients are linear in δk i,j and δl i,j . Unstable cases occur when N is singular, i.e. the determinant of N vanishes. In this case, small changes in δK induce large variations in δX. After expansion and factorisation, we obtain:

det(N) = 64(Z 2 Y 1 -Y 2 Z 1 )(Z 2 X 1 -X 2 Z 1 )(Y 2 X 1 -Y 1 X 2 ).
The unstable configurations may then be geometrically interpreted as the following two cases:

• the three camera planes are linearly dependent (i.e. the three projection planes in space do not intersect in a common point);

• all planes formed by any three space points intersect any image plane in a line parallel to the intersection lines of the camera planes.

This can be easily proved by assuming that i k is a unit vector otrthogonal to the plane P k with k ∈ {1, 2, 3}.

We first examine the case where the three vectors i k are linearly dependent. This is equivalent to the case of two images. The third image can be deduced from the two first ones.

The second case is when the three vectors i k are linearly independent. As the vectors {i 1 , i 2 , i 3 } can be defined as an orthogonal basis, it is clear that points M 0 , M 1 , M 2 have (0, 0, 0), (X 1 , Y 1 , Z 1 ), (X 2 , Y 2 , Z 2 ) as coordinates in this basis. The instability condition det(M ) = 0 is then equivalent to the fact that there exists k ∈ {1, . . . , 3} such that the projection of M 0 M 1 and M 0 M 2 on P k are collinear.

Structure and Motion Algorithm

The complete Structure and motion algorithm can now be described as follows. The key idea is to use polynomial resultants [START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF] to reduce the system to an univariate polynomial:

1. Computation of the coefficients. Compute the scale factors with the affine fundamental matrices [START_REF] Shapiro | 3D motion recovery via affine epipolar geometry[END_REF]. Obtain the coefficients k i,j and l i,j .

2. Elimination of the unknowns. Compute symbolically the resultant R 1 of P 1,2 and P 2,3 in Y 1 . Compute symbolically the resultant R 2 of R 1 and

P 1,3 in Y 2 . Set R 1,2 = √ R 2 (R 2 is a squared polynomial). Compute formally the resultant S 1 of Q 1,2 and Q 2,3 in Z 1 . Compute formally the resultant S 2 of S 1 and Q 1,3 in Z 2 . Set S 1,2 = √ S 2 (S 2 is a squared polynomial). Compute formally the resultant T of R 1,2 and S 1,2 in X 2 .
Finally, we obtain the polynomial T which is an univariate polynomial of degree 2 in X 2

1 . T = 0 can be easily solved in closed-form.

The algorithm gives directly the 3D coordinates of points in space, the projection matrices can then be computed from them.

Experimental results

Maximum Likelihood Estimation

The method described in this paper is minimal, in other words, there is no cost function and the results obtained are exact, depending only on input data. However, our algorithm involved solving a polynomial of degree 2 as well as a set of linear equations, which might induce numerical instabilities. For that purpose, we devise the Maximum Likelihood Estimator (MLE), for further refining the previously computed solution. In more detail, it consists in minimizing the Root Mean Squares of the reprojection residuals using an adequat parameterization of camera matrices and scene points. Let us denote the set of parameters as Θ. This set contains 12 motion parameters which are two rotations and translations and 3n structure parameters, where n is the number of points considered. Camera matrices are formed according to equations (1). The MLE is then given by argmin Θ C(Θ) where the cost function is defined as:

C(Θ) = i=3 i=1 j=n j=1 d 2 (m i j , mi j ),
where d 2 (., .) denotes the squared Metric distance in the image, m i j denotes points measured in the images and mi j denotes the reprojected points. The optimization is conducted using the Levenberg-Marquardt method [START_REF] Gill | Practical Optimization[END_REF].

We insist on the fact that, theoretically, both results should be the same. Comparing the solution provided by our algorithm and the Maximum Likelihood Estimator will only tells us about the numerical stability of the algorithm. The MLE is equivalent to bundle adjustment.

Simulated data

Minimal data. The simulation protocol is set up as follows. We choose at random 3 affine cameras, and four 3D points. We then compute the projections, and add noise on them. With the projections we can compute the 3D structure, which can be compared to the original points in space. This gives us one possible evaluation of reconstruction errors. This experiment is repeated 100 times with different noise levels from 0 to 10 pixels. Figure 2 shows the results. We can observe that the error is approximately linear in the added noise level.

RANSAC Estimation. We now take 157 3D points from a calibration object (Figure 4). We project them with computed projections, and then add some noise on the obtained images. Then, we compute a solution with RANSAC. First, we compare the results of the method with those obtained after a bundle adjustement on the motion computation. The bundle adjustment is made with the data computed by the method. Next, we do the same but the bundle adjustement is on structure and motion instead of only motion. The results illustrated in Figure 3 show that the obtained results are almost as good as the optimal solutions from bundle adjustements.

Real images

Target images. We carried out recontruction experiments from images of a target. With an appropriate software, we match points (there are no outliers in those images). The structure of the target gives us the affine calibration of the camera (We have seen that between the three images the calibration changes of less than 0.5%). The target is an object composed of three planes, two pairs of them forming a right angle (figure 4). At first we have seen that the average reprojection error is less than 0.5 pixel. In fact, because the reconstruction is Metric we expect to see the right angles between the planes in the reconstruction. Figure 5 shows the two top views of these angles.

House images. We have tested the algorithm on the house image sequence. We took three images of a house, detect points of interest with the Harris detector, match them with cross-correlation between adjacent images and triplet registration, see figure 6. We can see that there is a lot of bad matches. We first detect the outliers with a RANSAC method, so as to eliminate them. Next, we compute a model that minimizes the reprojection error over 100 trials. For the model, see figure 7.

Conclusion

We have presented a method for Metric reconstruction from calibrated affine cameras. The method works with minimal data, i.e. 4 points in 3 images. The stucture of 4 points in space is computed in closed-form rather than the motion parameters used by other researchers. We also studied the instability of minimal reconstruction. Experimental results show that the method is stable. Another important point is that the method is fast as it needs only to solve a polynomial of degree 2. As a comparison, bundle

Figure 2 .

 2 Figure 2. The median of the errors with different noise levels on simulated data.

Figure 3 .

 3 Figure 3. The two first graphs give a comparison between the new method followed by a bundle adjustement on motion parameters. The two next graphs give a comparison between the new method followed by a bundle adjustement on motion and structure. Note that the two first curves are undistinguishables.

Figure 4 .

 4 Figure 4. Three images of the target.

Figure 5 .

 5 Figure 5. The two right angles of the reconstructed target.

Figure 7 .

 7 Figure 7. The house with matched points.

Figure 6 .

 6 Figure 6. The house with matched points.