E Mouragnon

M Lhuillier
email: maxime.lhuillier@univ-bpclermont.fr

M Dhome

F Dekeyser

P Sayd

Generic and Real Time Structure from Motion using Local Bundle Adjustment 1

Keywords: Structure from Motion, Bundle Adjustment, Generic camera

This paper describes a method for estimating the motion of a calibrated camera and the three dimensional geometry of the filmed environment. The only data used is video input. Interest points are tracked and matched between frames at video rate. Robust estimates of the camera motion are computed in real-time, key frames are selected to enable 3D reconstruction of the features. We introduce a local bundle adjustment allowing 3D points and camera poses to be refined simultaneously through the sequence. This significantly reduces computational complexity when compared with global bundle adjustment. This method is applied initially to a perspective camera model, then extended to a generic camera model to describe most existing kinds of cameras. Experiments performed using real world data provide evaluations of the speed and robustness of the method. Results are compared to the ground truth measured with a differential GPS. The generalized method is also evaluated experimentally, using three types of calibrated cameras: stereo rig, perspective and catadioptric.

1 Introduction

Previous Work

Automatic estimation of 3D scene structure and camera motion from an image sequence (known as "Structure from Motion" or SfM) has been the subject of much investigation. Different camera models -pinhole, fish-eye, stereo, catadioptric, multicamera systems, etc -have been used in such studies. This is still a very active field of research, and several successful SfM systems currently exist [START_REF] Pollefeys | Automated reconstruction of 3D scenes from sequences of images[END_REF]1,[START_REF] Nister | Visual odometry[END_REF][START_REF] Lhuillier | A quasi-dense approach to surface reconstruction from uncalibrated images[END_REF][START_REF] Royer | Localization in urban environments: monocular vision compared to a differential gps sensor[END_REF].

Calculation Time vs. Accuracy A large number of dedicated algorithms (i.e. for given camera models) have been successfully developed and are now commonly used for perspective or stereo rig models [START_REF] Nister | Visual odometry[END_REF][START_REF] Pollefeys | Automated reconstruction of 3D scenes from sequences of images[END_REF]. There are several noteworthy types of dedicated SfM algorithms. These include firstly approaches without global optimization of the full geometry which are fast but of questionable accuracy (due to errors that accumulate over time). Among the work proposed for Vision-Based SLAM (Simultaneous Localization and Mapping), Nistér et al. [START_REF] Nister | Visual odometry[END_REF] have presented an approach known as "visual odometry". This method estimates the motion of a stereo head or a simple camera in real-time, using visual data only: its aim is to provide guidance for robots. Davison [3] proposes a real-time camera pose calculation method but assumes a small number of landmarks only (less than 100 landmarks). This approach is best suited to indoor environments and is not therefore appropriate for lengthy camera displacements due to the complexity of its algorithms and progressively increasing uncertainties. There are also different off-line methods that use bundle adjustment to optimize global geometry and thus obtain highly accurate models (see [START_REF] Triggs | Bundle adjustment -A modern synthesis[END_REF] for a thorough survey of bundle adjustment algorithms). Such an optimization is very costly in terms of computing time and can not be implemented in a realtime application. Bundle adjustment entails iterative adjustment for camera poses and point positions in order to obtain the optimal least squares solution. Most of these articles call on the Levenberg-Marquardt (LM) algorithm, which combines the Gauss-Newton algorithm with the method of gradient descent to solve the non linear criterion involved in bundle adjustment. The main problem in such bundle adjustment is that it is very slow, especially for long sequences because it requires inversion of linear systems whose size is proportional to the number of estimated parameters (even if the sparse structure of the systems involved is taken into account). It is also important to have an initial estimate relatively close to the real solution. Applying a bundle adjustment in a hierarchical way is therefore an interesting idea [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF][START_REF] Shum | Efficient bundle adjustment with virtual key frames: A hierarchical approach to multi-frame structure from motion[END_REF][START_REF] Royer | Localization in urban environments: monocular vision compared to a differential gps sensor[END_REF] but does not solve the computing time problem. An alternative method is therefore necessary to decrease the number of parame-ters to be optimized. Shum et al. [START_REF] Shum | Efficient bundle adjustment with virtual key frames: A hierarchical approach to multi-frame structure from motion[END_REF] exploit information redundancy in images by using two virtual key frames to represent a sequence. Steedly and Essa [START_REF] Steedly | Propagation of innovative information in Non-Linear Least-Squares structure from motion[END_REF] propose an incremental reconstruction with bundle adjustment. They then readjust only the parameters which have changed. While their method is faster than a global optimization, it is not efficient enough for long sequences that are highly data-dependent. Kalman filters or extended Kalman filters [START_REF] Davison | Real-time simultaneous localization and mapping with a single camera[END_REF] are another possibility, but such filters are known to provide less accurate results than bundle adjustment. Generic vs. Specific Structure from Motion Yet another widely explored avenue is that of omni-directional central (catadioptric, fish-eye) or non-central (e.g. multicamera) systems that offer a larger field of view [START_REF] Chang | Omni-directional structure from motion[END_REF][START_REF] Micusik | Autocalibration & 3D reconstruction with noncentral catadioptric cameras[END_REF][START_REF] Pless | Using many cameras as one[END_REF]. It is highly challenging to develop generic SfM tools suitable for any camera model. This approach has recently been investigated using generic camera models [START_REF] Kannala | A generic camera model and calibration method for conventional, wide-angle, and fish-eye lenses[END_REF][START_REF] Ramalingam | A generic structure-from-motion framework[END_REF]. In such models, pixels define image rays in the camera's coordinate system. Rays corresponding to pixels are given by the calibration function. They intersect at a single point usually called "projection center" in the case of a central camera [START_REF] Kannala | A generic camera model and calibration method for conventional, wide-angle, and fish-eye lenses[END_REF] but this is not necessary in other cases. In recent work on generic SfM, camera motion is estimated by the generalization of the conventional essential matrix [START_REF] Ramalingam | A generic structure-from-motion framework[END_REF][START_REF] Pless | Using many cameras as one[END_REF][START_REF] Mouragnon | Generic and real-time structure from motion[END_REF] derived from the Pless Equation [START_REF] Pless | Using many cameras as one[END_REF], and minimal relative pose estimation algorithms [START_REF] Stewénius | Solutions to minimal generalized relative pose problems[END_REF]. In the same way as for the specific models, a method is also required for refinement of 3D points and camera poses. In the generic case which implies different cameras (pinhole, stereo, fish-eye), bundle adjustment is different from the conventional approach used for perspective cameras. The minimized error may be a 3D or a 2D error. As the projection function is not explicit for some camera models, a 3D error can be used [START_REF] Ramalingam | A generic structure-from-motion framework[END_REF][START_REF] Micusik | Autocalibration & 3D reconstruction with noncentral catadioptric cameras[END_REF]. A 3D error is not however optimal since it favors the contribution of far points from the cameras and can produce biased results [START_REF] Lu | Fast and Globally Convergent Pose Estimation from Video Images[END_REF]. An other solution is to minimize a 2D reprojection error (in pixels) by clustering all camera rays such that each cluster of rays is approximated by a perspective camera [START_REF] Ramalingam | A generic structure-from-motion framework[END_REF].

Summary and comparison with previous work

This paper first proposes an accurate and fast incremental reconstruction and localization algorithm. Our idea [START_REF] Mouragnon | Real time localization and 3d reconstruction[END_REF] is to take advantage of both offline methods with bundle adjustment and faster incremental methods. In our algorithm, a local bundle adjustment is carried out each time a new camera pose is added to the system. A similar approach [START_REF] Engels | Bundle adjustment rules[END_REF] was also published a few months after ours [START_REF] Mouragnon | Real time localization and 3d reconstruction[END_REF], but it is not generic and it does not include a key-frame selection to stabilize 3D calculation as ours. A related approach is proposed by Zhang and Shan [START_REF] Zhang | Incremental motion estimation through modified bundle adjustment[END_REF], but their work calls for local optimization of an image triplet only, and eliminates structure parameters from the proposed reduced local bundle adjustment. By taking into account 2D reprojections of 3D estimated points in more than three images without eliminating 3D point parameters, it is possible to greatly improve reconstruction accuracy.

Secondly, this paper shows how our fast reconstruction method is extended to generic cameras. Within this generalization framework, our method [START_REF] Mouragnon | Generic and real-time structure from motion[END_REF] replaces image projections of a specific camera model with use of back-projected rays and minimization of angular error between rays. Its first advantage is, of course, a high degree of interchangeability between camera models. The second advantage is that it is also effective where the image projection function is not explicit (as in the case of non-central catadioptric cameras) and precludes clustering.

In summary, previously proposed approaches related to ours include the following:

generic but non-real-time methods [START_REF] Ramalingam | A generic structure-from-motion framework[END_REF][START_REF] Kannala | A generic camera model and calibration method for conventional, wide-angle, and fish-eye lenses[END_REF][START_REF] Lhuillier | Effective and Generic Structure from Motion[END_REF] ( [START_REF] Kannala | A generic camera model and calibration method for conventional, wide-angle, and fish-eye lenses[END_REF] deals only with central cameras).

real-time but non-generic methods [START_REF] Davison | Real-time simultaneous localization and mapping with a single camera[END_REF][START_REF] Nister | Visual odometry[END_REF], not using bundle adjustment. generic methods using the Pless Equation [START_REF] Pless | Using many cameras as one[END_REF][START_REF] Ramalingam | A generic structure-from-motion framework[END_REF] (generalization of the epipolar constraint), but without details on how to solve the equation in common situations.

dedicated or "specific" (non-generic) methods using local bundle adjustment [START_REF] Mouragnon | Real time localization and 3d reconstruction[END_REF][START_REF] Engels | Bundle adjustment rules[END_REF].

Our Contribution

The first new feature of our approach is use of local bundle adjustment, for both specific model (perspective) and generic model cameras. The second is inclusion of a detailed method for solving the Pless equation (in most cases, this is not a "simple" linear problem as suggested in [START_REF] Pless | Using many cameras as one[END_REF][START_REF] Ramalingam | A generic structure-from-motion framework[END_REF]). Our original study of artificial solutions of the linearized version of Pless equation is a byproduct. Finally, the system as a whole is new (in that it is the first to be both real-time and generic).

For purposes of clarity, this paper first describes the complete method (including local bundle adjustment) for a standard (perspective) camera in Section 2. This Section also compares the time complexity of our local bundle adjustment in the incremental scheme to that of the standard (global) bundle adjustment in the hierarchical scheme. Sections 3 and 4 deal with the generalization of the method, by describing the generic camera model and explaining modifications to geometry refinement. This is followed by a more detailed look at the generic initialization step and various solutions to the Pless equation. Finally, experiments performed for perspective and generic camera models are presented in Sections 5 and 6 respectively.

Incremental Method for the Perspective Camera Model

Let us consider a video sequence acquired with a camera mounted on a vehicle moving through an unknown environment. The purpose of our work is to enable determination of camera position and orientation in a global reference frame at several points in time t, along with the 3D position of a set of points (viewed along the scene). To do so, we use a monocular camera whose intrinsic parameters (including radial distortion) are known and assumed to be unchanged throughout the sequence. The algorithm initially determines a first image triplet for use in setting up the global frame and system geometry (Section 2.2). A robust pose calculation is then carried out for each frame of the video flow (Section 2.3) using feature detection and matching (Section 2.1). Some of the frames are selected as key frames for 3D point triangulation (Section 2.4). The system operates "incrementally", and when a new key frame and 3D points are added, it performs local bundle adjustment (Section 2.5). The result (see Figure 4) is a set of camera poses corresponding to key frames and 3D coordinates of the points seen in the images.

Interest Point Detection and Matching

Our entire method is based on the detection and matching of feature points (see Figure 1). In each frame, Harris corners [START_REF] Harris | A combined corner and edge detector[END_REF] are detected. A pair of frames is matched as follows:

• For each interest point in image 1, we select possible corresponding points in a region of interest defined in image 2. • Then a Zero Normalized Cross Correlation score is computed for these interest points neighborhoods. • The pairs with the highest scores are selected to provide a list of corresponding point pairs between the two images. A unicity constraint (winner takes all) is applied such that a point can only be matched with one other point.

For adaptation to a real-time application, the "detection and matching" step has been implemented using SIMD extensions of modern processors. This is both a fast and highly efficient solution.

Sequence Initialization With 3 Views

We know that the motion between two consecutive frames must be sufficiently large to compute the epipolar geometry and 3D points. For this reason, frames are selected relatively far from each other but with a suitable number of common points. To do so, the first image denoted I 1 is always selected as a key frame. The second image I 2 is selected as far as possible from I 1 in the video but with at least M interest points matched in I 1 . Then for I 3 , the frame selected is the one farthest from I 2 , so that there are at least M matched interest points between I 3 and I 2 and at least M ′ matched points between I 3 and I 1 (for our experiments, the values M = 400 and M ′ = 300 were selected). This process affords a sufficient number of corresponding points between frames to calculate the movement of the camera. The camera coordinate system associated with I 1 is taken as the global coordinate system, and the relative poses between the first three key frames are calculated using the 5-point algorithm [START_REF] Nister | An efficient solution to the five-point relative pose problem[END_REF] and a RANSAC [START_REF] Fischler | Random sample consensus, a paradigm for model fitting with applications to image analysis and automated cartography[END_REF] approach. Use of three views enhances the robustness of the method and eliminates ambiguities induced by coplanar points [START_REF] Nister | An efficient solution to the five-point relative pose problem[END_REF]. More details on the initialization process are given in [START_REF] Royer | Localization in urban environments: monocular vision compared to a differential gps sensor[END_REF]. Observed points are subsequently triangulated into 3D points using the first and the third observation. Finally estimated poses and 3D point coordinates are optimized through standard bundle adjustment

Robust Real-Time Pose Estimates

Let's assume that poses obtained with cameras C 1 to C i-1 and corresponding to selected key frames I 1 to I i-1 were previously calculated in the reconstruction reference frame. A set of points were computed whose projections were present in the corresponding images. The next step is to calculate camera pose C corresponding to the last acquired frame I. To do so, we match I (last acquired frame) and I i-1 (last selected key frame) to determine a set of points p whose projections on the cameras (C i-2 C i-1 C) are known and whose 3D co-ordinates have already been computed. Taking 3D points reconstructed from C i-2 and C i-1 , we use Grunert's pose estimation algorithm as described in [START_REF] Haralick | Review and analysis of solutions the three point perspective pose estimation problem[END_REF] to compute the location of camera C. A RANSAC process then gives an initial estimate of camera C pose which is subsequently refined using a fast LM optimization stage with only 6 parameters (3 for optical center position and 3 for orientation). At this stage, the covariance matrix of the camera pose C is calculated via the Hessian inverse and we can draw a 90% confidence ellipsoid (see Figure 2). If Cov is the covariance matrix of camera pose C, the ellipsoid of confidence is given by ∆x T Cov -1 ∆x ≤ 6.25 since ∆x T Cov -1 ∆x obeys the X 2 distribution with 3 DOFs.

Key Frame Selection and 3D Point Reconstruction

As already seen above, not all frames of the input are taken into account for 3D reconstruction: only a sub-sample of the video is used (Figure 3). For each frame, the normal approach is to compute the corresponding localization from the last two key frames. In our case, a criterion is incorporated to indicate whether or not a new frame must be added as a key frame. First, if the number of points matched with the last key frame I i-1 is not sufficient (typically, below a fixed level M, (M = 400 with in our experiments), a new key frame is required. This is also necessary if the calculated position uncertainty exceeds a certain level (for example, more than the mean distance between two consecutive, key positions). Obviously, the frame that did not meet the criterion cannot become a new key frame I i and the immediately preceding frame is therefore taken. New points (i.e. those only observed in I i-2 , I i-1 and I i ) are reconstructed using a standard triangulation method.

Local Bundle Adjustment

Following selection of the last key frame I i and its addition to the others, the reconstruction is optimized. The optimization process is a bundle adjustment or Levenberg-Marquardt minimization of the cost function f i (C i , P i ) where C i and P i are respectively the parameters of the cameras (extrinsic parameters) and the 3D points chosen for this stage i. The idea is to reduce the number of calculated parameters by optimizing only the extrinsic parameters of the n last cameras and accounting for the 2D reprojections in the N (with N ≥ n) last frames (see Figure 5). Thus, C i = {C i-n+1 . . . C i } and P i contains all the 3D points projected on cameras C i . The cost function f i is the sum of points P i reprojection errors in the last frames C i-N +1 to C i :

f i (C i , P i ) = C k ∈{C i-N+1 ... C i } P j ∈P i ||ε k j || 2 where ||ε k j || 2 = d 2 (p k j , K k P j )
is the square of the Euclidean distance between K k P j , estimated projection of 3D point P j through the camera C k and the measured corresponding observation p k j . K k is the projection matrix 3 × 4 of camera k comprising C k extrinsic parameters and known intrinsic parameters. Thus, n (number of optimized cameras at each stage) and N (number of images taken into account in the reprojection function) are the two main parameters involved in the optimization process. Their given value can influence both quality of results and speed of execution. Our experiments enabled the determination of values for n and N (typically, we use n = 3 and N = 10) that provide an accurate reconstruction. Optimization takes place in two LM stages with an update of inliers/outliers between the two stages. A series of LM iterations is stopped if the error is not suitably decreased or a maximum number of iterations is reached (5 in our case). In practice, the number of necessary iterations for each local bundle adjustment is quite low; this is due to the fact that, excepting the last added camera, all the cameras poses have already been optimized at stage i -1, i -2, ... It should be noted that when the reconstruction process starts, we refine not only the last parameters of the sequence, but the very whole 3D structure. Thus, for i ≤ N f , we opted for N = n = i. N f is the maximum number of cameras required for stage i optimization to be global (in our experiments, we choose N f = 20). In this way, reliable initial data is obtained, an important factor given the recursive nature of the algorithm, without any problems, since the number of parameters is still relatively restricted at this stage.

Time Complexities for Local and Global Bundle Adjustments

Detailed calculation complexity calculations are given in the Appendix A of this paper. Let p be the number (considered as constant) of points projected through each camera. Based on a video sequence of N seq key frames, the complexity of one local bundle adjustment iteration applied to the n last cameras after allowance for 2D reprojections on N images is

C i C i-1 C i-2 C i-3 n N C i P i
Θ(p.N + p.n 2 + n 3 )
whereas the complexity of one global bundle adjustment iteration applied to the whole sequence is Θ p.N 2 seq + N 3 seq .

It is clearly advantageous to reduce the number of parameters (n and N) involved in the optimization process. As an example, the complexity gains with respect to global bundle adjustment obtained for a sequence of 20 key frames and 150 2D reprojections per image are shown in Table 3. This study deals with only one LM iteration. In practice, parameters n and N are fixed with low values (n = 3 and N ≤ 10) for our incremental method. Furthermore, the number of iterations is limited by 10 for each local bundle adjustment (although we note that the mean number of necessary iterations is less than that). As a consequence, the time complexity of local bundle adjustment is O(1) for each key frame and is O(N seq ) in our incremental scheme for a complete sequence.

This should be compared to the complexity of global bundle adjustment in the hierarchical scheme [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. At first glance, one iteration of the last bundle adjustment step over the whole sequence is greater than N 3 seq . So it is easy to say that the time complexity is the best for the incremental scheme. A more valuable comparison could be done if we assume that the track length is bounded by a constant l. In that case, the complexity N 3 seq due to solve the reduced camera system of global bundle is reduced to N seq l 2 [START_REF] Triggs | Bundle adjustment -A modern synthesis[END_REF]. So the complexity of one iteration of global bundle is at least proportional to the sequence length. Since the global bundle adjustment is applied in a hierarchical scheme, we apply it on 1 sequence of length N seq , 2 sequences of length N seq /2, 4 sequences of length N seq /4 ... (in the reverse order). Thus the overall complexity is greater than N seq log(N seq ). One time again, the time complexity is the best for the incremental scheme.

Method Summary

In summary, the proposed method consists of the following steps:

(1) Select an image triplet that provides the first three key frames (Section 2.2). Set up the global frame, estimate the relative pose, and triangulate 3D points. (2) For each new frame, calculate matches with last key frame (Section 2.1) and estimate the camera pose and uncertainty (Section 2.3). Determine whether a new key frame is needed. If not, repeat 2. (3) If a new key frame is necessary, select the preceding frame as new key frame, triangulate new points (Section 2.4) and make a local bundle adjustment (Section 2.5). Repeat the above starting from step 2.

Generic Camera Model and Geometry Refinement

The method presented in Section 2 is designed for a perspective camera model. Our approach, however, is geared to also using other kinds of cameras for 3D reconstruction. Dedicated methods are possible for each (e.g. catadioptric camera, stereo rig) but a method suitable for any type of camera involved would be very valuable. Our approach thus entails extending the incremental method to a generic camera model [START_REF] Mouragnon | Generic and real-time structure from motion[END_REF]. In the following section (Section 3), the pixel error applied to refine geometry is replaced by a generic error. Robust initialization with three views is then described in detail in Section 4.

Camera Model

For any pixel p of a generic image, the (known) calibration function f of the camera defines an optical ray r = f (p). This projection ray is an oriented line r = (s, d) for which s is the starting point or origin and d is the direction of the ray in the camera frame (||d|| = 1). For a central camera, s is a single point (camera center) whatever the pixel p. In the general case, s could be any point given by the calibration function.

Error Choice

Let P j = [x j , y j , z j , t j ] ⊤ be the homogeneous coordinates of the j-th point in the world frame. Let R i and t i be the orientation (rotation matrix) and the origin of the i-th camera frame in the world frame.

If (s i j , d i j ) is the optical ray corresponding to the observation of P j through the i-th camera, the direction of the line defined by s i j and P j is

D i j = R i ⊤ [I 3 | -t i ]P j -t j s i
j in the i-th camera frame. In the ideal case, directions d i j and D i j are parallel (which is equivalent to an image reprojection error of zero pixels).

The conventional approach [START_REF] Triggs | Bundle adjustment -A modern synthesis[END_REF][START_REF] Mouragnon | Real time localization and 3d reconstruction[END_REF] consists of minimizing a sum of squares ||ǫ i j || 2 where ǫ i j is a specific error depending on the camera model: the reprojection error in pixels. In our case, a generic error must be minimized. We thus define ǫ i j as the angle between the directions d i j and D i j defined above (see Figure 6).

Many experiments show that the convergence of bundle adjustment is very poor with ǫ i j = arccos(d i j .

D i j ||D i j ||
) and satisfactory with ǫ i j defined as follows [START_REF] Lhuillier | Effective and Generic Structure from Motion[END_REF]. We have thus chosen

ǫ i j = π(R i j D i j ) with R i j a rotation matrix such that R i j d i j = [0 0 1] ⊤ and π a function R 3 → R 2 such that π([x y z] ⊤ ) = [ x z y z ] ⊤ . Note that ǫ i
j is a 2D vector whose Euclidean norm ||ǫ i j || is equal to the tangent of the angle between d i j and D i j . The tangent is a good approximation of the angle if it is small.

Intersection, Resection and Bundle Adjustment

Once ǫ i j redefined, it is a straightforward matter to redefine the common tools for (incremental) reconstruction. Calculation of a 3D point coordinates P j knowing its observations in cameras C 1 . . . C n (intersection) is given by min-Camera i(R i , t i ) P oint j(x j , y j , z j , t j ) D i j d i j ǫ i j s i j Fig. 6. Angular bundle adjustment: the angle between observation ray (s i j , d i j ) and 3D ray D i j which goes from s i j to 3D point is minimized.

imizing g j :

g j (P j ) = i=1...n ||ǫ i j || 2 .
Calculation of camera pose C i knowing 3D points P 1 . . . P m (resection) is given by minimizing g i :

g i (C i ) = j=1...m ||ǫ i j || 2 .
In the same way as for the "specific" method, a local bundle adjustment is applied each time a new key frame I i is added to the reconstruction. Parameters of estimated 3D points and cameras at the end of the sequence are refined by the minimization of the cost function g i :

g i (C i , P i ) = C k ∈{C i-N+1 ... C i } P j ∈P i ||ǫ k j || 2 .
where C i and P i are respectively the generic camera parameters (extrinsic parameters of key frames) and the 3D points chosen for stage i (Figure 7). We account for points reprojections in the N (with N ≥ n) last frames (typically n = 3 and N = 10).

Generic Initialization

The following paragraphs describe in detail the generic initialization of incremental 3D reconstruction. Section 4.1 briefly presents the Plücker coordinates used to describe 3D rays in space. The generic epipolar constraint (or Pless Equation) is also described in Section 4.2, and Section 4.3 shows how to solve it in different cases. Finally robust initialization with three views and robust pose estimation are explained in Sections 4.4 and 4.5 respectively.

C i C i-1 C i-2 C i-3 P i C i N n
Fig. 7. Local angular bundle adjustment when camera C i is added. Only points P i and cameras C i parameters surrounded by dotted lines are optimized. The minimized criterion nevertheless accounts for 3D points projections in the N last images.

Plücker Coordinates

Based on a set of pixel correspondences between two images, the relative pose (R, t) of two cameras can be estimated in a generic framework. For each 2D points correspondence (x 0 , y 0 ) and (x 1 , y 1 ) between images 0 and 1, we have a correspondence of optical rays (s 0 , d 0 ) and (s 1 , d 1 ). A ray (s, d) is defined by its Plücker coordinates, which are convenient for this calculation. The Plücker coordinates of a 3D line (or ray) are two 3 × 1 vectors (q, q ′ ) respectively named direction vector and moment vector. q gives the direction of the line and q ′ is such that q • q ′ = 0. Any point P on the line described by (q, q ′ ) satisfies q ′ = q ∧ P. Plücker coordinates are defined up to a scale and a parametrization of the line is given by (q ′ ∧ q) + αq, ∀α ∈ R if ||q|| = 1.

For an optical ray (s, d) as previously defined, where s is the origin of the ray and d its direction in the camera frame, Plücker coordinates in the same frame are:

     q = d q ′ = d ∧ s

Generalized Epipolar Constraint (or Pless Equation)

Let camera 0 be the origin of the global coordinates system. Its pose is written as R 0 = I 3 and t 0 = 0 0 0 ⊤ . We want to determine (R, t) which is the pose of camera 1 in the same frame. For each pixel correspondence between image 0 and image 1, optical ray (q 0 , q ′ 0 ) and (q 1 , q ′ 1 ) must intersect in space at a single 3D point M (see Figure 8).

Camera 0 = global f rame M (q 0 , q ′ 0 ) Camera 1 (R, t) (q 1 , q ′ 1 )
Fig. 8. Relative pose of 2 cameras

In the global frame, the two rays are:

     q 00 = q 0 q ′ 00 = q ′ 0 and      q 10 = Rq 1 q ′ 10 = -[t] × Rq 1 + Rq ′ 1 where [t] × =        0 -t z t y t z 0 -t x -t y t x 0        and t = t x t y t z ⊤ .
[t] × is the skew-symmetric cross-product matrix of t such that [t] × x = t ∧ x for a 3 × 1 vector x.

The two rays intersect if and only if q 10 • q ′ 00 + q ′ 10 • q 00 = 0, which leads to the generalized epipolar constraint or Pless Equation [START_REF] Pless | Using many cameras as one[END_REF]:

q ′ 0 ⊤ Rq 1 -q ⊤ 0 [t] × Rq 1 + q ⊤ 0 Rq ′ 1 = 0 (1)
The generalized essential matrix G is defined as follows [START_REF] Ramalingam | A generic structure-from-motion framework[END_REF]:

G =    -[t] × R R R 0 3×3   
The matrix G is a 6 × 6 matrix verifying the constraint:

   q 1 q ′ 1    ⊤ G    q 0 q ′ 0    =
0 if and only if rays (q 0 , q ′ 0 ) and (q 1 , q ′ 1 ) intersect in space. G contains 9 zero coefficients and the 9 coefficients of R twice. There are thus 18 useful coefficients.

We identify two cases where this equation has an infinite number of solutions. Obviously, this number is infinite if the camera is central (the 3D is recovered up to a scale). We note that Equation 1 is the usual epipolar constraint defined by the essential matrix E = [t] × R if the camera center is at the origin of the camera frame.

The second case is less obvious but it occurs in practice. In our experiments, we assume there are only simple matches: all projection rays (s i , d i ) of a given 3D point go through a same camera center (in the local coordinate of the generic camera). In other words, we have q ′ 0 = q 0 ∧ c 0 and q ′ 1 = q 1 ∧ c 1 with c 0 = c 1 . For a multicamera system comprising central cameras (e.g. the stereo rig), this means that 2D points correspondences are only made with points of the same sub-image. This is often the case in practice due to the small regions of interest used for reliable matching, or the empty intersections between the fields of views of compositing cameras. If the camera motion is a pure translation (R = I 3 ), Equation 1 becomes q ⊤ 0 [t] × q 1 = q ′ 0 ⊤ q 1 + q ⊤ 0 q ′ 1 = 0 where the unknown is t. In this context, the scale of t can not be estimated. We assume for our purposes that the camera motion is not a pure translation at the initialization stage.

Solving the Pless Equation

Equation 1 is rewritten as

q ′ 0 ⊤ Rq 1 -q ⊤ 0 Ẽq 1 + q ⊤ 0 Rq ′ 1 = 0 (2)
where the two 3 × 3 matrices ( R, Ẽ) are the new unknowns. We store the coefficients of ( R, Ẽ) in an 18 × 1 vector x and see that each value of the 4-tuple (q 0 , q ′ 0 , q 1 , q ′ 1 ) produces a linear equation a ⊤ x = 0. If we have 17 different values of this 4-tuple for each correspondence k, we have 17 equations a ⊤ k x = 0. This is enough to determine x up to a scale factor [START_REF] Ramalingam | A generic structure-from-motion framework[END_REF]. We have built the matrix A 17 containing the 17 correspondences such that

A 17 x = 0 with A ⊤ 17 = [a ⊤ 1 |a ⊤ 2 | • • • a ⊤ 17 ]
. The resolution depends on the dimension of the A 17 kernel which directly depends on the type of camera used. We determine Ker(A 17 ) and its dimension by a Singular Value Decomposition of A 17 . In this paper, we have distinguished three cases: (1) central cameras with an unique optical center (2) axial cameras with collinear centers and (3) nonaxial cameras.

It is not surprising that the kernel dimension of the linear system to solve is greater than one. Indeed, the linear Equation 2 has more unknowns (18 unknowns) than the non-linear Equation 1 (6 unknowns). Possible dimensions are reported in Table 2 and are justified below. Previous works [START_REF] Pless | Using many cameras as one[END_REF][START_REF] Ramalingam | A generic structure-from-motion framework[END_REF] ignored these dimensions, although a (linear) method is heavily dependent on them. Central Camera For central cameras (e.g. pinhole cameras), all optical rays converge at the optical center c. Since

q ′ i = q i ∧ c = [-c] × q i , Equation 2 becomes q 0 ⊤ ([c] × R -Ẽ -R[c] × )q 1 = 0. We note that ( R, Ẽ) = ( R, [c] × R -R[c] ×
) is a possible solution of equation 2 for any 3 × 3 matrix R. Such solutions are "exact": Equation 2 is exactly equal to 0 whatever (q 0 , q 1 ). Our "real" solution is ( R, Ẽ) = (0, [t] × R) if c = 0, and it is not exact due to image noise. Thus the dimension of Ker(A 17 ) is at least 9 + 1. Experiments have confirmed that this dimension is 10 (up to noise). In this case, we simply solve the usual epipolar constraint constraint q 0 ⊤ [t] × Rq 1 = 0 as described in [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF].

Axial Camera This case includes the common stereo rig of two perspective cameras. Let c a and c b be two different centers of the camera axis. Appendix B shows that "exact" solutions ( R, Ẽ) are defined by

Ẽ = [c a ] × R -R[c a ] × and R ∈ V ect{I 3×3 , [c a -c b ] × , (c a -c b )(c a -c b ) ⊤ }
based on our assumption of "simple" matches (Section 4.2). Our real solution is not exact due to image noise, and we note that the dimension of Ker(A 17 ) is at least 3+1. Experiments have confirmed that this dimension is 4.

We build a basis of 3 exact solutions x 1 , x 2 , x 3 and a non-exact solution y with the singular vectors corresponding to the four smallest singular values of A 17 . The singular values of x 1 , x 2 , x 3 are 0 (up to computer accuracy) and that of y is 0 (up to image noise). We calculate the real solution ( R, Ẽ) by linear combination of y, x 1 , x 2 and x 3 such that the resulting matrix R verifies R⊤ R = λI 3×3 or Ẽ is an essential matrix. Let l be the vector such that l ⊤ = [λ 1 λ 2 λ 3 ] ⊤ , and thus we denote as R(l) and Ẽ(l) the matrix R and Ẽ extracted from solution

y -[x 1 |x 2 |x 3 ]l. Using these notations, we have R(l) = R 0 -3 i=1 λ i R i and Ẽ(l) = E 0 -3 i=1 λ i E i with (R i , E i ) extracted from x i .
Once the basis x 1 , x 2 , x 3 is calculated, we compute the coordinates of the solution by non-linear minimization of the function (λ, l) → λI 3×3 -R(l) ⊤ .R(l) 2 to obtain l and thus Ẽ. SVD decomposition is applied to Ẽ, and we obtain 4 solutions [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] for ([t] × , R). The solution with the minimal epipolar constraint A 17 x is then selected. Lastly, we refine the 3D scale k by minimizing

k → i (q ′ 0i ⊤ Rq 1i -q ⊤ 0i k.[t] × Rq 1i + q ⊤ 0i Rq ′ 1i ) 2 and perform t ← kt.
Non-Axial Camera For a non-axial camera (e.g. a multicamera system with perspective cameras such that centers are not collinear), the problem is also different. Appendix B shows that the "exact" solutions are ( R, Ẽ) ∈ V ect{(I 3×3 , 0 3×3 )} based on our assumption of "simple" matches (Section 4.2). Our real solution is not exact due to image noise, and we see that the dimension of Ker(A 17 ) is at least 1+1. Experiments have confirmed that this dimension is 2. We have not yet experimented this case.

Initialization with Three Views (RANSAC process)

The first step of the incremental algorithm is the 3D reconstruction of a subsequence containing the first key frames triplet {0, 1, 2}. A number of random samples are taken, each containing 17 points. For each sample, the relative pose between views 0 and 2 is computed using the above described method and matched points are triangulated. The pose of camera 1 is estimated with 3D/2D correspondences by iterative refinement minimizing the angular error (see details in Section 3.3). The same error is minimized to triangulate points. Finally, the solution producing the highest number of inliers in views 0, 1, and 2 is selected from among all samples. The j-th 3D point is considered as an inlier in view i if the angular error ||ǫ i j || is less than ǫ (ǫ = 0.01 rad in our experiments).

Pose Estimates (RANSAC)

The generic pose calculation is useful for both steps of our approach (initialization and incremental process). We assume that the i-th pose P i = (R i , t i ) of the camera is close to that of the i-1-th pose P i-1 = (R i-1 , t i-1 ). P i is estimated by iterative non-linear optimization initialized at P i-1 with a reduced sample of five 3D/2D correspondences, in conjunction with RANSAC. For each sample, the pose is estimated by minimizing an angular error (Section 3.3) and we count the number of inliers (points) for this pose. The pose with the maximum number of inliers is then selected and another optimization is applied with all inliers.

Experiments using a Perspective Camera

We have applied our incremental localization and mapping algorithm to a semi-urban scene. The goal of these experiments was to evaluate robustness (resistance to perturbations) in a complex environment and compare accuracy to the ground truth provided by a Real-Time Kinematics Differential GPS. The camera was therefore mounted on an experimental vehicle whose velocity was about 1.1 m/s. The distance covered was about 70 m and the video sequence was 1 min long (Image size was 512 × 384 pixels at 7.5 f ps). More than 4.000 3D points were reconstructed and 94 images selected as key frames from a series of 445 (figure 11). This sequence was particularly interesting because of image content (Figure 9 shows people walking in front of the camera, sunshine, etc...) which was not conducive to the reconstruction process. Moreover, the environment was more appropriate to a GPS localization (which provides ground truth) since the satellites involved were not occluded by high buildings.

Processing Time

Calculations were performed on a standard Linux PC (Pentium 4 at 2.8 GHz, 1Go of RAM memory at 800 MHZ). Image processing times through the sequence were as shown in Figure 10. Times measured included feature detection (#1500 Harris points per frame), matching, and pose calculation for all frames. For key frames, longer processing times were necessary (see Figure 10) due to 3D points reconstruction and local bundle adjustment. We took n = 3 (number of optimized camera poses) and used N = 10 (number of cameras for minimization of reprojection criterion). Note that computing speeds were interesting with an average of 0.09 s for normal frames and 0.28 s for key frames. These data are shown in Table 3.

Frames

Comparison with Ground Truth (Differential GPS)

The calculated trajectory obtained with our algorithm was compared against data given by a GPS sensor. For purposes of comparison, we applied a rigid transformation (rotation, translation and scale factor) to the trajectory as described in [START_REF] Faugeras | The representation, recognition, and locating of 3-d objects[END_REF] to fit with GPS reference data. Figure 12 shows trajectory registration with the GPS reference. As GPS positions are given in a metric frame we could compare camera locations and measure positioning error in meters. For camera key pose i, 3D position error is:

E i3D = (x i -x GP S ) 2 + (y i -y GP S ) 2 + (z i -z GP S ) 2
and 2D position error in the horizontal plane is:

E i2D = (x i -x GP S ) 2 + (y i -y GP S ) 2
where x i , y i , z i are the estimated coordinates for camera pose i and x GP S , y GP S , z GP S are the corresponding GPS coordinates. Figure 13 shows 2D/3D error variations for the 94 frames. The maximum measured error was 2.0 m with a 3D mean error of 41 cm and a 2D mean error of less than 35 cm.

Parametrizing Local Bundle Adjustment

In our incremental method, local bundle adjustment consists of optimizing the end of the reconstruction only, to avoid unnecessary calculations and excessive computing times. As mentioned before, optimization is applied to the n last estimated camera poses, with allowance for points reprojections in a larger number N of cameras. We therefore tested several values for n and N as reported in Tables 4, 5 and 6.

First, note that we must have N ≥ n+2 to define the reconstruction frame and the scale factor at the end of the sequence. For N = n or n + 1, it happened that the reconstruction was not completed due to process failure before the end of the sequence (Tables 4 and5). This confirms that N ≥ n + 2 is an important condition. We also measured mean time processing for local bundle adjustment as a function of n (Table 6). In practice, it does not vary much with N since the mean track length of points in key-frames (Figure 14) is limited.

Then, we compare results with the GPS trajectory (Table 4) and a trajectory computed using global bundle adjustment (Table 5). As expected, the quality increases when n increase. We also observed that higher values of n provides minor improvements of quality in our context. In all the following experiments, we fix n = 3 and N = 10 since this provides a good trade-off between accuracy and time performance.

Comparison with Global Bundle Adjustment for a Long Sequence

An other experiment has been carried out in a town-center (Figure 15) where differential GPS is not available for ground truth. One can visually ensure that reconstruction is not much deformed and drift is low compared to the covered distance. The incremental method selects 354 key-frames among 2900 frames and reconstructs 16135 3D points with n = 3 and N = 10. The estimated mean 3D position error compared to global bundle adjustment is 29cm (the trajectory length is about 500m). That shows that our algorithm (very appropriate to long scene reconstruction in term of computing time) is also robust and gives result similar as that of global bundle adjustment.

Experiments using a Generic Camera

The incremental generic 3D reconstruction method was tested on real data with three different cameras: a perspective camera, a catadioptric camera and a stereo rig. Examples of frames are shown in Figure 16 and sequences characteristics in Table 7. Computing performances are reported in Table 8. Following these experiments, the trajectory obtained with our generic method was compared to GPS ground truth or global specific bundle adjustment results. As already seen above, a rigid transformation (rotation, translation and scale factor) must be applied to the trajectory to fit with reference data [START_REF] Faugeras | The representation, recognition, and locating of 3-d objects[END_REF]. A mean 3D error or 2D error in the horizontal plane can then be measured between the generic and the reference trajectory. The same matching method was used for all three.

Comparison with Ground Truth (Differential GPS)

The first results correspond to the same video sequence as in Section 5.2, with a pinhole camera mounted on an experimental vehicle equipped with a differential GPS receiver (inch precision). The calculated motion obtained with our generic algorithm was compared against data given by the GPS sensor, and Figure 17 shows the two recorded trajectories. As GPS positions are given in a metric frame, we could compare camera locations and measure positioning error in meters: mean 3D error was 0.48 m and 2D error in the horizontal plane was 0.47 m. Results were slightly less accurate than those obtained with the specific perspective method.

Comparison with Specific and Global Bundle Adjustment

In the following examples, ground truth is not available. We therefore compare our results against those of the best method available: a global and specific bundle adjustment (all 3D parameters have been refined so as to obtain an optimal solution with a minimum reprojection error). Sequences characteristics and results are reported in Table 7.

Sequence 1 is taken in an indoor environment using an handheld pinhole camera. A very small difference is obtained: the mean 3D error is less than 6.5 cm for a trajectory length of (about) 15 m. The relative error is 0.45%. Sequence 2 is taken in an outdoor environment, using an handheld catadioptric camera (the 0-360 mirror with the Sony HDR-HC1E camera shown on Figure 18, DV format). The useful part of the rectified image is contained in a circle whose diameter is 458 pixels. The difference is also small: the mean 3D error is less than 9.2 cm for a trajectory length of (about) 40 m. Relative error is 0.23%. Sequence 3 is taken with a stereo rig (baseline: 40 cm) in a corridor (Figure 18). The image is composed of two sub-images of 640 × 480 pix. The trajectory (20 m long) is compared to results obtained using the left/right camera and global bundle adjustment. Mean 3D error is 2.7/8.4 cm compared with the left/right camera and relative error is 0.13/0.42%.

Conclusion

We have developed and tested a method for solving the real-time Structure from Motion problem. This paper presents a complete process for doing so, from feature detection and matching to estimating geometry and refining it by local bundle adjustment. 3D reconstruction is an incremental process that begins with initialization using three key frames. For subsequent frames, camera localization is then determined and key frames are selected for 3D point reconstruction. This method is extended to a generic camera model facilitating the change from one kind of camera to another. Promising results have been obtained on real data with three different types of cameras (even if there is room for further optimizing implementation). Results and time complexities are compared favorably with those of global bundle adjustment (and ground truth when available). Future work includes experiments on more complex multicamera systems and automatic 3D modeling methods using the generic camera model.

Appendix A: Complexity of One Bundle Adjustment Iteration

Let P be the vector of parameters to be estimated (orientation of the cameras + position of their optical centers + 3D points coordinates), X the set of 3D points projections detected in images and f (P) the projection of 3D points in images according to the parameters we are looking for. The problem here is to minimize the function φ(P) = ||f (P) -X|| 2 . At stage k of the iterative algorithm, we calculate ∆ k such that P k+1 = P k + ∆ k . ∆ k is obtained by solving the equation: J ⊤ J.∆ k = J ⊤ .ǫ k where J is the Jacobian matrix of f calculated in P k and ǫ k = X -f (P k ) (more precisely, the diagonal terms of matrix J ⊤ J are multiplied by a coefficient using the Levenberg-Marquardt method [START_REF] Press | Numerical recipes in C[END_REF]).

Since we want to process long sequences with multiple parameters to be evaluated (six for each camera and three for each 3D point), it seems natural to make use of bundle adjustment characteristics, i.e. the block structure of matrix J ⊤ J for reconstruction of a set of points [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. This matrix is composed of three blocks U, V, and W such that U and V are block-diagonal:

• U, matrix made of diagonal 6 × 6 blocks representative of the dependency of image i measurements on the associated camera parameters. • V, matrix made of diagonal 3 × 3 blocks representative of relations between point j parameters and the measurements associated with this point. • W, matrix expressing the inter-correlations between 3D points parameters and cameras parameters. The structure of W is linked to the fact that many

V j i W ⊤
U W Fig. 19. Structure of the J ⊤ J matrix points are not projected through all cameras. W has a number of non-null 6 × 3 blocks equal to the number of 2D reprojections (half of the dimension of X or f (P)).

The system is therefore

   U W W ⊤ V       ∆ cameras ∆ points    =    Y cameras Y points  
, and it is solved in two steps [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]:

(1) Calculation of the increment ∆ cameras to be applied to cameras by resolution of the following system:

(U -WV -1 W ⊤ )∆ cameras = Y cameras -WV -1 Y points (3)
(2) Direct calculation of the increment ∆ points applicable to 3D points:

∆ points = V -1 (Y points -W ⊤ ∆ cameras )
Let C and P be the number of cameras and points optimized in bundle adjustment. Let p be the number (considered as constant) of points projected through each camera.

Once J ⊤ J is calculated (Figure 19) (time complexity is proportional to the number N r of 2D reprojections taken into account), the two most time-consuming expensive stages of this resolution are:

• calculation of matrix product WV -1 W ⊤ • resolution of camera linear system (3).

For matrix product WV -1 W ⊤ , the number of necessary operations can be determined by first considering the number of non-null blocks of WV -1 . This is the same number as W, i.e. (p.C), number of reprojections in C images, because V -1 is block-diagonal. Then, in the product (WV -1 )W ⊤ , each non-null 6×3 block of WV -1 is used once in the calculation of each block column of WV -1 W ⊤ . Thus the time complexity of the product WV -1 W ⊤ is Θ(p.C 2 ). The time complexity

Fig. 1 .

 1 Fig. 1. Example of an image from video data. Small squares represent detected interest points, and white lines the distance covered by the matched points.

Fig. 2 .

 2 Fig. 2. Top view of a processing reconstruction. This shows the trajectory, reconstructed 3D points and a confidence ellipsoid for the most recently calculated camera.

Fig. 3 .

 3 Fig. 3. Video sub-sampling: localization of all frames. 3D point reconstruction with key frames (squares).

Fig. 4 .

 4 Fig. 4. Top view of a complete reconstruction in an urban environment. The distance covered is a length of about 200 meters including a half-turn. More than 8.000 3D points are reconstructed for 240 key frames.

Fig. 5 .

 5 Fig. 5. Local bundle adjustment when camera C i is added. Only the points and cameras surrounded by dotted lines are optimized. We nevertheless account for 3D point reprojections in the N last images.

Fig. 9 .

 9 Fig. 9. Two frames used for the real data experiments.

Fig. 10 .

 10 Fig. 10. Processing time (in seconds) for non-key frames and key frames.

Fig. 11 .

 11 Fig. 11. Top view of the reconstructed scene and trajectory (# 4.000 points and 94 key positions).

Fig. 12 .

 12 Fig. 12. Registration with the GPS reference, top: in horizontal plane, bottom: along altitude axis. The continuous line represents GPS trajectory and points the estimated key positions. Coordinates are expressed in meters.

Fig. 13 .

 13 Fig. 13. Error in meters. Continuous line = 2D error, dotted line = 3D error.

Fig. 14 .

 14 Fig. 14. Point distribution with track length.

Fig. 15 .

 15 Fig. 15. Two images of the urban video, map and top view of the reconstructed scene and trajectory (16135 points and 354 key positions).

Fig. 16 .

 16 Fig.[START_REF] Mouragnon | Generic and real-time structure from motion[END_REF]. Feature tracks for a generic camera image, using three types of cameras: perspective (top left), catadioptric (top right), and stereo rig (bottom) cameras. The same matching method was used for all three.

Fig. 17

 17 Fig. 17. Top: Registration of generic vision trajectory with GPS ground truth. The continuous line represents GPS and points represent vision estimated positions. Bottom: 3D error along the trajectory.

  Fig. 17. Top: Registration of generic vision trajectory with GPS ground truth. The continuous line represents GPS and points represent vision estimated positions. Bottom: 3D error along the trajectory.

Fig. 18 .

 18 Fig. 18. Cameras and top views of 3D reconstructions. Top: catadioptric camera (Sequence 2). Bottom: stereo rig (Sequence 3). Trajectory in blue and 3D points in black.

Table 1

 1 Complexity gain obtained through local bundle adjustment in comparison with global bundle adjustment for one iteration.

	Type	p	n N gain
	global	150 20 20	1
	reduction 1 150 5 20 10
	reduction 2 150 3 10 25

Table 2

 2 dim(Ker(A 17 )) depends on the kind of camera.

	Camera	Central Axial Non-Axial
	dim(Ker(A 17 ))	10	4	2

Table 6

 6 Mean 3D position error (in meters) compared to GPS for the incremental method with different n and N values, and for global bundle adjustment. Mean local bundle adjustment computing times as a function of n for many N values (in seconds).

	H H H H H H n N n	n+1	n+2 n+3 n+5 n+7
	n=2	failed failed 0.55 0.49 0.85 1.99
	n=3	failed 3.28	0.45 0.41 0.41 0.41
	n=4	6.53	1.77	0.42 0.40 0.41 0.27
	global			0.33
	Table 4
	H H H H H H n N n	n+1	n+2 n+3 n+5 n+7
	n=2	failed failed 3.17 0.43 1.56 1.80
	n=3	failed 3.61	1.60 0.43 0.30 0.47
	n=4	7.67	1.44	1.03 0.24 0.25 0.36
	Table 5
	Mean 3D position error (in meters) compared to global bundle adjustment for
	different n and N values.

Table 7

 7 Characteristics of video sequences.

	Sequence		Camera	#Frames	#Key frames #3D Pts	#2D Pts	Traj. length
	Sequence 1		pinhole		511	48	3162	11966	15 m
	Sequence 2	catadioptric		1493	132	4752	18198	40 m
	Sequence 3	stereo rig		303	28	3642	14189	20 m
	Camera		Image size	Detection+Matching	Frame	Key frame Mean rate
	Pinhole		512 × 384		0.10	0.14	0.37	6.3 fps
	Catadioptric	464 × 464		0.12	0.15	0.37	5.9 fps
	Stereo rig		1280 × 480		0.18	0.25	0.91	3.3 fps

Table 8

 8 Computing times in seconds for the three cameras.

of the traditional resolution of the linear system (1) is Θ(C 3 ) [START_REF] Press | Numerical recipes in C[END_REF]. Therefore the time complexity of one bundle adjustment iteration is Θ(N r + p.C 2 + C 3 ).

Appendix B: Exact Solutions of the (linearized) Pless Equation

The i-th pair of rays is defined by Plucker Coordinates (q 0i , q ′ 0i ) and (q 1i , q ′ 1i ) such that q

The rays of the i-th pair intersect if they satisfy the Pless Equation

where the two 3 × 3 matrices ( R, Ẽ) are the unknowns. In this section, we seek an ( R, Ẽ) such that ∀i,

In other words, we seek an R ∈ ℜ 3×3 such that this Ẽ is independent of any available camera center pair (c 0i , c 1i ). As a consequence, Equation 4 is exactly equal to 0 (up to computer accuracy) whatever (q 0i , q 1i ). We consider many cases.

Simple Matches Only: ∀i,

As mentioned in Section 4.2, this particular case is important in practice.

Central Camera Let c be the center. This case is straightforward: we have

Stereo Camera Let c a and c b be the centers. We have

Furthermore, it is easy to show that this constraint does not allow another R by changing the coordinate basis such that c a -c b = 0 0 1

Axial Camera All camera centers of an axial camera are collinear: there are c a and c b such that

This case is the same as the previous one.

Non-Axial Camera There are three non-collinear centers c a , c b , c c . The constraint on R is

This constraint is three times that obtained for the stereo camera. Let's change the coordinate basis such that {c a -c b , c b -c c , c c -c a } is the canonical basis and write the solutions for the three stereo cases: any R ∈ V ect{I 3×3 } is possible.

Whether Matches are Simple or Not

In this case, the "simple match" constraint is not enforced.

Central Camera It does not occur here.

This constraint is stronger than that for the stereo camera with simple matches: we see that any

Axial Camera There are c a and c b such that c 0i = (1 -λ 0i )c a + λ 0i c b and c 1i = (1 -λ 1i )c a + λ 1i c b with λ 0i ∈ ℜ and λ 1i ∈ ℜ. Thus,

should not depend on λ 0 and λ 1 : we have 0

This case is the same as the previous one.

Non-Axial Camera There are three non-collinear centers c a , c b , c c and the constraint on R is three times the one obtained for the stereo camera:

Let's change the coordinate basis such that {c a -c b , c b -c c , c c -c a } is the canonical basis and write the solutions for the three stereo cases: we obtain R = 0. This method does not provide exact solutions.