Alexandre Eudes

Maxime Lhuillier

Error Propagations for Local Bundle Adjustment

Local bundle adjustment (LBA) has recently been introduced to estimate the geometry of image sequences taken by a calibrated camera. Its advantage over standard (global) bundle adjustment is a great reduction of computational complexity, which allows real-time performances with a similar accuracy. However, no confidence measure on the LBA result such as uncertainty or covariance has yet been introduced. This paper introduces statistical models and estimation methods for uncertainty with two desirable properties: (1) uncertainty propagation along the sequence and (2) real-time calculation. We also explain why this problem is more complicated than it may appear at first glance, and we provide results on video sequences.

Introduction

The robust and automatic estimation of camera motion and scene points from an image sequence (Structure-from-Motion or SfM) is still today an active field of research. A real-time method based on local bundle adjustment (LBA) was introduced [START_REF] Mouragnon | Real time localization and 3d reconstruction[END_REF] for this problem three years ago. However, no confidence measure on the result such as uncertainty or covariance was introduced before. A such measure would be useful to give a quality information if ground truth is not available, or even to merge the reconstruction with data provided by other sensors such as GPS or odometer. This Section summarizes previous works and our contribution on this topic.

Bundle Adjustment Bundle adjustment is a well known iterative method [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF] designed to solve non-linear least square problems for SfM. It estimates a vector x which minimizes a cost function x → ||y -F (x)|| 2 . Vector x usually concatenates 6 parameters for each camera pose and 3 parameters for each 3D point. Vector y concatenates detected features in images, F concatenates projection functions, and ||.|| is the Euclidean norm. We call this method global bundle adjustment (GBA) if x contains all parameters of the sequence.

The time complexity of one GBA iteration is O(c 3 + cp) with c and p the numbers of camera poses and 2D points, respectively. Although this complexity benefits the sparse structure of the problem and the assumption c ≪ p, realtime GBA is impossible for long sequences.

Error Propagation Error propagation provides confidence measures for bundle adjustment result [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. Assume that the image vector y follows a known Gaussian noise. Then, the function which maps y to minimizer x is approximated by its linear Taylor expansion. Now, the Gaussian noise y is propagated to Gaussian noise x such that the x covariance matrix can be estimated. Thanks to this matrix, uncertainty ellipsoids for a given probability can be defined for the bundle adjustment result.

The time complexity to estimate covariances of all camera poses and 3D points is at least that of one GBA iteration.

Extended Kalman Filter Bundle adjustment is based on the Levenberg-Marquardt method [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], which is an improved (damped) version of the Gauss-Newton Method. On the other hand, the Extended Kalman Filter (EKF) may be viewed as an incremental version of the Gauss-Newton Method [START_REF] Bertsekas | Incremental least squares methods and the extended kalman filter[END_REF] and has been applied for real-time SfM [START_REF] Davison | Real-time simultaneous localization and mapping with a single camera[END_REF]. EKF is known to provide less accurate results than bundle adjustment, but it is faster. It has is own covariance management for the estimated geometry (geometry and its covariance are jointly updated over time). EKF requires prior knowledge on the covariance of all estimated parameters. A recent improvement [START_REF] Bedder | Incremental estimation without specifying a-priori covariance matrices for the novel parameters[END_REF] cancels this drawback, but it still depends on linearization of constraints between frames.

Local Bundle Adjustment

We recently introduced LBA for real-time SfM [START_REF] Mouragnon | Real time localization and 3d reconstruction[END_REF]. The incremental SfM method is summarized as follows. Camera poses and a sparse cloud of 3D points are reconstructed for video frames before time t, and a new video frame should be added to the reconstruction at time t. Interest points [START_REF] Harris | A combined corner and edge detector[END_REF][START_REF] Lowe | Distinctive image features from scales-invariant keypoints[END_REF] of frame t are detected and matched with those of time t -1. Since many points of frame t -1 are reconstructed, there are 2D-3D correspondences for points in frame t and the new pose can be estimated. This pose is initialized by a robust method (such as RANSAC [START_REF] Fischler | Random sample consensus: a paradigm for model fitting with application to image analysis and automated cartography[END_REF]), and LBA is applied to refine the geometry (camera poses and 3D points) of the n most recent frames tn + 1, • • • , t -1, t by maintaining consistency with previous frames t-N +1, • • • , t-n (n < N ). LBA is a bundle adjustment such that [START_REF] Bay | Surf: Speeded up robust features[END_REF] x concatenates all 3D parameters of the n most recent frames and (2) y contains detected features in the N last frames.

The time complexity of one LBA iteration is O(n 3 + np) with p the numbers of 2D points in the N last images. Thanks to fixed and small values of n and N (e.g. n = 3, N = 10), the GBA complexity is greatly reduced.

Our Contribution

We introduce statistical models and covariance estimation methods for the SfM parameters estimated by LBA.

Section 3 describes a first method derived from the original LBA definition using a strong independence hypothesis (between the poses of previous frames t -N + 1, • • • , tn and the detected features in the N last frames). However, the original LBA is not a Maximum Likelihood Estimator (MLE) for this hypothesis. Then Section 4 describes a second method derived from a new LBA which is MLE. These two methods satisfy both desirable properties (1) uncertainty propagation along the sequence and (2) real-time performance.

Section 5 introduces covariance methods for both original and new LBAs using a weak independence hypothesis, but the real-time performance is lost. Last, we provide covariance results on video sequences and compare them with the covariance derived from GBA (Section 6).

Notations and Definitions

Time and Frame Index The integer t is both a time index and a frame index. The LBA-based SfM method distinguishes key-frames and non key-frames of the given video sequence to ensure a stable estimation of 3D. Since LBA is only applied on key-frames, we ignore non key-frame in this paper and use t as a key-frame number and time.

Camera Poses As mentioned above, LBA is bundle adjustment such that estimated parameters x are all 3D parameters of the n most recent frames and data y contains detected points in the N most recent frames. We assume that n = 3 and N = 10 to simplify notations without loss of generality. Notice that the camera pose of a frame is estimated many times by many LBAs since the sliding window size n is greater than one. For this reason, a double index is used for camera poses. After the LBA of time t, the poses of frames 0, 1

T ] T the vector which concatenates x and y. We define

x t = [c t-2 t c t-1 t c t t | • • • s i t • • • ] (1)
y t = [• • • m i t ′ • • • ] (2)
p t = [c t-9 t c t-8 t • • • c t-3 t ] (3)
such that LBA is concisely written as Other Notations Here we present our last definitions. A sub-vector of vector x is a vector obtained by removing one or many coordinates of x. A sub-matrix of a square matrix C is a matrix obtained from C by removing a line and a column of the same index, and this operation may be repeated for many indexes. We introduce

x t = f t (
x c t = [c t-2 t c t-1 t c t t ] as a sub- vector of x t . Notation z ∼ N (z, C z ) means that z is a Gaussian vector of mean z and covariance C z . Let z → g(z)
be a C 1 continuous function with Jacobian ∂g ∂z . Up to the first order, we have error propagation [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]

g(z) ∼ N (g(z), C g ) with C g = ∂g ∂z (z)C z ( ∂g ∂z (z)) T . ( 7
)
Recurrence Relation Now, the LBA-based SfM method may be written as a recurrence relation on camera poses. We estimate

[p t x c t ] from [p t-1 x c t-1 ] as follows: 1. p t is a sub-vector of [p t-1 x c t-1 ]
2. estimate x t = f t (p t , y t ) using LBA (Eq. 5)

3.

x c t is a sub-vector of x t Step 1 is obvious since [p t-1 x c t-1 ] = [c t-10 t-1 • • • c t-1 t-1 ] and p t = [c t-9 t • • • c t-3 t ] = [c t-9 t-1 • • • c t-3 t-1 ].

Original LBA and its Covariance

The goal of this Section is to link the covariance of [p t x c t ] to the recurrence relation of Section 2: we assume that

[p t-1 x c t-1 ] ∼ N ([p t-1 xc t-1 ], C [pt-1x c t-1 ] ) with a known covariance C [pt-1x c t-1 ]
and we would like to estimate

C [ptx c t ] such that [p t x c t ] ∼ N ([p t xc t ], C [ptx c t ]
). We assume that y t ∼ N (ȳ t , σ 2 I) with noise scale σ > 0 estimated by GBA during the SfM initialization [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF].

At first glance, the covariance C xt of the parameters x t estimated by LBA may be approximated by the inverse of approximated Hessian of LBA cost function at x t [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. However, this estimation does not propagate noise from the previous parameters p t to the new parameters x t ; it only propagates y t noise to x t noise. The practical consequence is that the uncertainty of camera poses will not grow with time, although this is the expected result. Our propagation methods do not have this problem, but they are more complicated.

Statistical Model of

[p t y t ] Covariance C pt is a sub-matrix of C [pt-1x c t-1 ] since p t is a sub-vector of [p t-1 x c t-1 ]
. Furthermore, we have y t ∼ N (ȳ t , σ 2 I). We assume that Gaussian vectors p t and y t are independent and obtain

[p t y t ] ∼ N ([p t ȳt ], C [ptyt] ), C [ptyt] = C pt 0 0 σ 2 I . (8)

First-Order Error Propagation

We approximate f t by its linear Taylor expansion at point [p t ȳt ] and obtain

p f (p, y) = p f (p, ȳ) + I 0 ∂f ∂p ∂f ∂y p - p y -ȳ . ( 9
)
Index t is omitted in this expression. We deduce that

[p t x t ] ∼ N ([p t xt ], C [ptxt] ) with covariance C [ptxt] = I 0 ∂ft ∂pt ∂ft ∂yt C [ptyt] I ∂ft ∂pt T 0 ∂ft ∂yt T . (10)
T ∂F ∂x ) -1 ∂F ∂x T (11) ∂f ∂p ≈ -( ∂F ∂x T ∂F ∂x ) -1 ∂F ∂x T ∂F ∂p (12)
with derivatives of F taken at point (p, f (p, y)).

Proof is provided in Appendix A (index t is omitted). Note that Eq. 12 is a generalization of the well known Eq. 11. The target covariance matrix

C [ptx c t ] is a top-left sub-matrix of C [ptxt] .

Algorithm

The LBA recurrence relation is completed with covariance. We estimate

[p t x c t ] and C [ptx c t ] from [p t-1 x c t-1 ] and C [pt-1x c t-1 ] as follows: 1. p t is a sub-vector of [p t-1 x c t-1 ] 2. C pt is a sub-matrix of C [pt-1x c t-1 ] 3. estimate x t = f t (p t ,
C [ptxt] = C pt C pt ∂ft ∂pt T ∂ft ∂pt C pt ∂ft ∂pt C pt ∂ft ∂pt T + σ 2 ∂ft ∂yt ∂ft ∂yt T (13)
The derivatives of f t are given by Eq. is the top-left sub-matrix of C [ptxt] such that

C [ptx c t ] = C pt A T A B , B is a top-left sub-matrix of σ 2 H -1 t + H -1 t ∂F t ∂x t T ∂F t ∂p t C pt ∂F t ∂p t T ∂F t ∂x t H -1 t and A is a top block of -H -1 t ∂F t ∂x t T ∂F t ∂p t C pt .
We estimate A and B with block-wise inversion [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]

H t = U W W T V ⇒ H -1 t = 0 0 0 V -1 + Y T Z -1 Y (14) with Z = U -WV -1 W T and Y = I -WV -1 . (15)
Matrices B, U and Z have the same dimensions; V is a block diagonal matrix with 3 × 3 invertible blocks thanks to Eq. 6. We successively estimate Z -1 , D 0 = ∂Ft ∂xt T ∂Ft ∂pt , and

D 1 = H -1 t D 0 using D 1 = 0 0 0 V -1 D 0 + (Y T (Z -1 (YD 0 ))). ( 16
)
At this point, the estimation of A and B is straightforward

: matrix A is a top block of -D 1 C pt and B-σ 2 Z -1 is a top-left sub-matrix of D 1 C pt D T 1 .

New LBA and its Covariance

In Section 3, the covariance for original LBA [START_REF] Mouragnon | Real time localization and 3d reconstruction[END_REF] is estimated assuming that p t and y t are independent. Section 4 introduces a new LBA function f n t which provides the maximum likelihood estimation (MLE) of x t for this assumption. Section 4 also describes a recurrence relation to estimate geometry and its covariance:

[p t x c t ] and C [ptx c t ] are estimated from [p t-1 x c t-1 ] and C [pt-1x c t-1 ] .

Maximum Likelihood Estimator

Index t is omitted for p t , pt , x t , xt , y t , F t , f t in this part. Assume that p ∼ N (p, C p ) and y ∼ N (ȳ, σ 2 I) are independent with ȳ = F (p, x). The unknown parameters of the statistical model are p and x (we assume that σ 2 and C p are given as true values). The probability density function of Gaussian vector [py] is

d(p, y|p, x) = Ke -1 2 ( 1 σ 2 ||y-F (p,x)|| 2 +(p-p) T C -1 p (p-p)) (17)
with K a constant. Thus, the MLE

f n t (p, y) of [px] is argmin [px] 1 σ 2 ||y -F (p, x)|| 2 + (p -p) T C -1 p (p -p). (18)
Function f n t defines estimations p n t and x n t of pt and xt (note that f t does not provide an estimation of pt ).

First-Order Error Propagation

Function f n t (p t , y t ) in Eq. 18 may be rewritten as

f n t (y n t ) = argmin xn ||y n t -F n t (x n t )|| 2 (19) with xn t = [p t |x t ], y n t = [(C pt ) -1 2 p t | 1 σ y t ] (20)
and

F n t (x n t ) = [(C pt ) -1 2 pt | 1 σ F t (p t , xt )]. (21)
Now we see that the Jacobian ∂f n t ∂y n t of f n t can be estimated using Eq. 11 of Proposition 1:

∂f n t ∂y n t = ( ∂F n t ∂ xn t T ∂F n t ∂ xn t ) -1 ∂F n t ∂ xn t T . (22)
Furthermore, Eq. 20 implies y n t ∼ N (ȳ n t , I). Eq. 22 and the linear Taylor expansion of f n t at [p t ȳt ] provides

C [p n t x n t ] = C f n t =
∂F n ∂x n = ∂F n ∂p ∂F n ∂x = C -1 2 p 0 1 σ ∂F ∂p 1 σ ∂F ∂x (24)
C [p n t x n t ] = C -1 p + 1 σ 2 ∂F ∂p T ∂F ∂p 1 σ 2 ∂F ∂p T ∂F ∂x 1 σ 2 ∂F ∂x T ∂F ∂p 1 σ 2 ∂F ∂x T ∂F ∂x -1 . ( 25
)
We apply block-wise inversion Eqs. 14 and 15 to C [p n t x n t ]

with x = x c x s , U = C -1 p + 1 σ 2 ∂F ∂p T ∂F ∂p 1 σ 2 ∂F ∂p T ∂F ∂x c 1 σ 2 ∂F ∂x c T ∂F ∂p 1 σ 2 ∂F ∂x c T ∂F ∂x c (26)
W = 1 σ 2 ∂F ∂p ∂F ∂x c T ∂F ∂x s , V = 1 σ 2 ∂F ∂x s T ∂F ∂x s (27)
and obtain C [ptx c t ] = (U -WV -1 W T ) -1 . Thank to Eq. 6, V is a 3 × 3-block diagonal matrix and is easy to inverse.

Weak Hypothesis

Previous Sections 3 and 4 require that p t ∼ N (p t , C pt ) and y t ∼ N (ȳ t , σ 2 I) are independent. However, camera pose c t-3 t = c t-3 t-1 is a sub-vector of both p t and x t-1 such that x t-1 = f t-1 (p t-1 , y t-1 ). Since y t and y t-1 have common 2D points in frames t -9, t -8, • • • t -1, we can not assert that p t and y t are independent.

In Section 5,

Statistical model of [y t p t ]

New notations are needed here. Let y t ∩ y t-1 (respectively, y t \ y t-1 ) be the sub-vector of y t with 2D points which are (respectively, which are not) in y t-1 . Covariance C [yt-1pt-1x c t-1 ] is known, and covariance C [yt∩yt-1|pt] is a sub-matrix of C [yt-1pt-1x c t-1 ] since y t ∩y t-1 is a sub-vector of y t-1 and p t is a sub-vector of [p t-1 x c t-1 ]. We assume that y t \ y t-1 and [y t ∩ y t-1 |p t ] are independent. Thus [y t p t ] is Gaussian vector with covariance

C [ytpt] = C [yt\yt-1|yt∩yt-1|pt] = σ 2 I 0 0 C [yt∩yt-1|pt] . (28)

First-Order Error Propagation (Original LBA)

We approximate f t by its linear Taylor expansion at point [ȳ t pt ] and obtain

  y p f (p, y)   =   ȳ p f (p, ȳ)   +   I 0 0 I ∂f ∂y ∂f ∂p   y - ȳ p -p . ( 29
)
Index t is omitted in this expression. We deduce that

[y t p t x t ] ∼ N ([ȳ t pt xt ], C [ytptxt]
) with covariance

C [ytptxt] =   I 0 0 I ∂ft ∂yt ∂ft ∂pt   C [ytpt] I 0 ∂ft ∂yt T 0 I ∂ft ∂pt T . (30)
Derivatives of f t are provided by Proposition 1.

Algorithm (Original LBA)

Sections 5.2 and 5.1 define the recurrence relation for the original LBA. We estimate [y t p t x c t ] and

C [ytptx c t ] from [y t-1 p t-1 x c t-1 ] and C [yt-1pt-1x c t-1 ] as follows: 1. p t is a sub-vector of [y t-1 p t-1 x c t-1 ] 2. C [yt∩yt-1|pt] is a sub-matrix of C [yt-1pt-1x c t-1 ]
3. estimate x t = f t (p t , y t ) using LBA (Eq. 5)

d(p, y|p, x) = Ke -1 2 zT C -1 [yp] z, z = y -F (p, x) p -p (31)
with K a constant. Thus, the MLE

f n t (p, y) of [px] is argmin [px] zT C -1 [yp] z with z = y -F (p, x) p -p . ( 32
)
Function f n t defines estimations p n t and x n t of pt and xt (note that Eq. 18 is a special case of Eq. 32).

First-Order Error Propagation (New LBA)

Algorithm (New LBA)

Sections 5.4 and 5.5 define the recurrence relation for the new LBA (Eq. 32). We estimate [y t p t x c t ] and

C [ytptx c t ] from [y t-1 p t-1 x c t-1 ] and C [yt-1pt-1x c t-1 ] as follows: 1. p t is a sub-vector of [y t-1 p t-1 x c t-1 ] 2. C [yt∩yt-1|pt] is a sub-matrix of C [yt-1pt-1x c t-1 ]
3. estimate C [ytpt] using Eq. 28 4. estimate [p n t x n t ] = f n t (p t , y t ) using LBA (Eq. 32) 5. estimate C [ytp n t x n t ] using Eqs. 37, 11 and 35 6. do p t ← p n t , x t ← x n t and

C [ytptxt] ← C [ytp n t x n t ] 7. C [ytptx c t ] is a top-left sub-matrix of C [ytptxt] .
6. Experiments

Integrating Covariance to LBA-based SfM

Our real-time SfM system [START_REF] Mouragnon | Real time localization and 3d reconstruction[END_REF] has two steps: initialization and incremental reconstruction. The former estimates the camera poses and 3D points of the sequence beginning using standard GBA. The latter incrementally reconstructs the sequence (poses and points) using original LBA (Eq. 5).

Then we integrate our covariance methods in the incremental step. These methods also require covariance for the camera poses at the sequence beginning. This covariance is estimated by the standard method derived from GBA [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]: the inverse of approximated Hessian of the minimized cost function, multiplied by image noise σ 2 .

A simple gauge is chosen to estimate the covariance of poses at the beginning: we fix the first frame pose (R 0 , t 0 ) with rotation R 0 = I and location t 0 = 0 and the largest coordinate of the t 0 -th frame location t t0 with t z t0 = 1 (t 0 = 9). This information should be given since it is known that the shape of uncertainty ellipsoids derived from covariance highly depends on the gauge choice [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]. Then we remove the columns and rows corresponding to these 7 parameters in the approximated Hessian before inversion.

How to Check LBA-based Covariance ?

LBA-based SfM method produces geometry estimations which are similar to those of GBA-based SfM [START_REF] Mouragnon | Real time localization and 3d reconstruction[END_REF]. So we expect to obtain the same result for geometry covariance: LBA-based covariance (our methods) should be similar to GBA-based covariance (standard method). We will compare both.

The last step of "GBA-based" SfM is GBA for the complete sequence: the vector x of all 3D parameters minimizes the cost function x → ||y -F (x)|| 2 , where y is the vector of all tracked 2D points along the whole sequence and F concatenates the corresponding projection functions. If y ∼ N (ȳ, σ 2 I), we have the GBA-based covariance

C x = σ 2 ( ∂F ∂x T ∂F ∂x ) -1 [8]
. This is the inverse of the approximated Hessian with the same gauge choice as the GBA for sequence beginning (Section 6.1). The value of σ is also the same. The complete calculation of Hessian inverse is not necessary: we only calculate the diagonal blocks we need thanks to Eqs. 14 and 15 (U and V are the sub-hessians of camera poses and 3D points [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF][START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF], respectively).

Results

Figure 1 shows three images of the sequence taken in urban area. The camera is calibrated, is mounted on a car and is pointing forward. The trajectory length is about 400m and the sequence has 2731 512 × 384 images. 384 key-frames are selected from the video. 16365 points are reconstructed by original LBA from 74236 Har- ris points [START_REF] Harris | A combined corner and edge detector[END_REF] matched using SURF descriptor [START_REF] Bay | Surf: Speeded up robust features[END_REF] in keyframes (we slightly modify the original SURF method and re-implement it on GPU using CUDA for real-time performance). The means of points in a frame and track lengths are 193 and 4.5, respectively.

Figure 2 shows quantitative comparisons of GBA-based covariance with LBA-based covariance described in Section 5.3 (original LBA with weak independence hypothesis). The same comparison is made in Figure 3 for the LBA-based covariance described in Section 4 (new LBA with strong independence hypothesis). In both cases, we study the major axis of the uncertainty ellipsoid of the location of camera c t-2 t with probability 90% (the t-2 th camera is updated at times t-2,t-1,t due to our sliding windows size of LBA, and we choose the uncertainty at the last update). The x-axis is the key-frame number; the sequence beginning optimized by GBA is not considered in these figures.

The tops of Figures 2 and3 show the ratio of major axis lengths between LBA and GBA. We see that the ratio is acceptable (close to 1.1) for the original LBA. Unfortunately, this method is not real-time. Furthermore, the ratio of new LBA is small (close to 0.6) due to the strong independence hypothesis between p t and y t . The bottoms of Figures 2 and3 show the angle between LBA and GBA major axes. The angles are acceptable (small) for both LBAs.

We have also experimented the original LBA-based covariance with strong hypothesis (Section 3). In this case, the ratio of major axis lengths diverges. Therefore the strong hypothesis should not be used with the original LBA.

At this point, the covariance of new LBA with strong hypothesis is the only choice in our real-time context (although its scale is too small). Let ē and σ e be the mean and standard deviation of ratios of major axis lengths between original LBA (weak hypothesis) and new LBA (strong hypothesis) for all key-frames of the sequence. We estimate ē = 1.82 and σ e = 0.13. Since σ e /ē is low, we decide to improve new LBA covariances by multiplying them with ē2 . Now, the main axis lengths of new LBA ellipsoids are roughly the same as those of original LBA. Figure 4 shows a top view of the reconstructed sequence with uncertainty ellipsoids of camera locations of our amended covariances. We see that 1. ellipsoids shapes are similar for new LBA and GBA 2. the major axis length increases progressively with time These are expected results. Only 6.4 ms are needed by the new LBA-based covariance method for each key-frame, so our method is real-time. The total time of new LBA covariance is 2.4 s, which is (obviously) smaller than that of GBA covariance (145 s) for all camera poses.

Our experiments also includes Monte-Carlo simulations (to check our implementation of LBA covariances) and covariance estimations with other real sequences (similar results are obtained with similar values of ē and σ e ).

Conclusion

This paper has introduced four covariance estimation methods for Structure-from-Motion (SfM) based on local bundle adjustment (LBA). They are derived from two noise hypotheses ("weak" and "strong") and two LBAs: our new LBA which provides Maximum Likelihood Estimation for these hypotheses, and the original LBA which does not. All methods propagate uncertainty along the sequence. Only two of them are real-time thanks to the strong hypothesis.

We must find a pair (statistical model, estimator) such that the estimated covariance is both physically plausible and real-time. On one side, realistic statistical models may require too much calculation to obtain real-time performance. On the other side, unrealistic statistical models may produces real-time but unrealistic uncertainty estimation. We have experimented our covariance methods on real sequences reconstructed by LBA-based SfM and have compared the results with standard covariance of global bundle adjustment. The original LBA with weak hypothesis provides acceptable covariance. This hypothesis is realistic but it does not allow real-time performance. Furthermore, the original LBA can not be used with the strong hypothesis. The new LBA with strong hypothesis provides acceptable covariance if empirical coefficient is introduced. This last method is the only choice in our real-time context.

Future work includes error propagation from key-frames to non key-frames, integration of calibration uncertainty in error propagation, experiments for many gauge choices, and fusion of our vision results with GPS or odometer data.

Appendix A

Proposition 1 (Section 3.2) is a particular case of Proposition 6.1 in [START_REF] Faugeras | Geometry of Multiple Images[END_REF], which also asserts that function f locally exists and is C 1 continuous. In this appendix, the time index t is omitted. Furthermore, F k , x i and p j are the k-th, i-th and j-th coordinates of vectors F , x and p, respectively.

We drop y for Eq. 12 proof since y acts a constant embedded in F . Second order partial derivatives of function

We obtain Eq. 12 since ∂F ∂x has full rank. Eq. 11 is obtained here in the special case F (p, x) = F (x)p by swapping notations y and p.

Vector p t fixes the gauge (coordinate frame and scale) for LBA since p t concatenates at least two camera poses [START_REF] Mouragnon | Real time localization and 3d reconstruction[END_REF]. In this context, ∂F ∂x has full rank for general configurations of 3D points [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF].

Acknoledgements This paper was supported by CNRS, CEA, ANR and Num@tec Automotive within the framework of ODIAAC project. We thank E. Mouragnon, S. Lion and S. Naudet for many technical improvements of the structure-from-motion software, and G. Jacob for GPU implementations of SURF and Harris using CUDA.

Figure 1 .

 1 Figure 1. Three images of the sequence.

Figure 2 .

 2 Figure 2. Top: ratio of major axis lengths between original LBA (weak hypothesis) and GBA ellipsoids. Bottom: angle between major axes of original LBA and GBA ellipsoids.

Figure 3 .

 3 Figure 3. Top: ratio of major axis lengths between new LBA (strong hypothesis) and GBA ellipsoids. Bottom: angle between major axes of new LBA and GBA ellipsoids.

Figure 4 .

 4 Figure 4. Top view of the reconstructed sequence by the LBAbased SfM. Camera locations and 3D reconstructed points are black dots. The 90% uncertainty ellipsoids of camera locations are also drawn for key-frames whose numbers are multiple of 5 (continuous line for standard GBA and dotted line for our new LBA with strong independence hypothesis).

  The first step is initialization and fixes the coordinate frame of the whole reconstruction. The second step is incremental and does not modify the pose of frames t-9, t-8, • • • t-3. In other words, ∀t ′ ≤ t -3, c t ′ t = c t ′ t-1 .

	SfM method has two steps:
	-∃t 0 > 0 such that c 0 t0 , c 1 t0 , • • • c t0 t0 are estimated by GBA.
	-∀t > t 0 , c t-2 t c t-9 t-1 , c t-8 t-1 , • • • c t-3 , c t-1 t t-1 .	, c t t are estimated by LBA from
	Local Bundle Adjustment LBA of time t also estimates
	the set of 3D points {s i t } which have at least one detected projection in frames t -2, t -1 and t. Let m i t ′ be the detec-tion of point s i t in frame t ′ such that t ′ ∈ {t -9, t -8, • • • t} (if any). Vectors c t ′ t , s i t , m i t ′ have dimensions 6,3 and 2, re-
	spectively. If x and y are vectors, we define [xy] = [x|y] =
	[x T
	, • • • t are defined by vectors c 0 t , c 1 t , • • • c t t respec-tively. Vector c t ′ t does not exist if t < t ′ . The LBA-based

  p t , y t ) with (4) f t (p t , y t ) = argmin xt ||y t -F t (p t , xt )|| 2 . (5)

	Function F t concatenates projections functions of 3D points
	{s i t } onto 2D points {m i t ′ } up to image noise. Different
	notations are used here for minimizer x t and variable vector xt of function F t . Thus, f t minimizes [c t-2 t ct-1 t ct t |∀s i t ] →
		||m i t ′ -c t ′ t	si t || 2 +	||m i t ′ -ct ′ t	si t || 2 (6)
	∀i,9≤t-t ′ ≤3		∀i,2≤t-t ′ ≤0
	with c t ′ t	si t , ct ′ t	si t projections of si t by poses c t ′ t and ct ′ t .

  y t ) using LBA (Eq. 5) 4. estimate C [ptxt] using Eqs. 8, 10, 11, and 12 5. C [ptx c t ] is a top-left sub-matrix of C [ptxt] . We only estimate block C [ptx c t ] of C [ptxt] since the dimension of camera vector [p t x c t ] is quite smaller than that of [p t x t ] which also contains 3D points parameters. Now, technical details are given for this estimation. Combining Eq. 10 and 8 we obtain

  We may estimate [p t x c t ] and C [ptx c t ] from [p t-1 x c t-1 ] and C [pt-1x c t-1 ] as follows: 1. p t is a sub-vector of [p t-1 x c t-1 ] 2. C pt is a sub-matrix of C [pt-1x c ] is a top-left sub-matrix of C [ptxt] . The new LBA (Eq. 18) is more time consuming than the original LBA (Eq. 5) since camera poses of p n t should be estimated with those of x n t . So we replace step 3 by x t = f t (p t , y t ) and remove step 5. Step 4 is unchanged: Eqs. 23 and 21 are still used with p n t = p t , x n t = x t , C [ptxt] = C [p n t x n t ] . In other words, we approximate the result of the new LBA by the result of the original LBA (both results are exactly the same if there is no image noise). We only estimate block C [ptx c t ] of C [ptxt] since the dimension of camera vector [p t x c t ] is quite smaller than that of [p t x t ] which also contains 3D points parameters. Now, technical details are given for this estimation. Eqs. 21 and 23 provide (t and ˜omitted)

								4.3. Algorithm
								t-1 ]
								3. estimate [p n t x n t ] = f n t (p t , y t ) using LBA (Eq. 18)
								4. estimate C [p n t x n t ] using Eqs. 23 and 21
								5. do p t ← p n t , x t ← x n t and C [ptxt] ← C [p n t x n t ]
								6. C [ptx c t
	∂f n t ∂y n t	I	∂f n t ∂y n t	T	= (	∂F n t ∂ xn t	T ∂F n t t ∂ xn	) -1 . (23)

  recurrence relations are introduced to estimate covariances of original and new LBAs without this hypothesis: [y t p t x c t ] and C [ytptx c t ] are estimated from [y t-1 p t-1 x c t-1 ] and C [yt-1pt-1x c t-1 ] .

. Maximum Likelihood Estimator (New LBA)

  Index t is omitted for p t , pt , x t , xt , y t , F t , f t in this part. Assume that [yp] ∼ N ([ȳ p], C [yp] ) with C [yp] defined by Eq. 28 and ȳ = F (p, x). The unknown parameters of the statistical model are p and x (we assume that C [yp] is given as true value). The probability density function of the Gaussian vector [py] is

	5.4
	4. estimate C [ytptxt] using Eqs. 28, 30, 11 and 12
	5. C [ytptx c t ] is a top-left sub-matrix of C [ytptxt] .
	Now the full estimate of ∂ft ∂yt by Eq. 11 is required for the
	top-left sub-matrix C [ytptx c t ] of C [ytptxt] and the method is
	intractable for real-time.


Maxime Lhuillier

Toward Automatic 3D Modeling of Scenes using a Generic Camera Model

The automatic reconstruction of 3D models from image sequences is still a very active field of research. All existing methods are designed for a given camera model, and a new (and ambitious) challenge is 3D modeling with a method which is exploitable for any kind of camera. A similar approach was recently suggested for structure-frommotion thanks to the use of generic camera models. In this paper, we first introduce geometric tools designed for 3D scene modeling with a generic camera model. Then, these tools are used to solve many issues: matching errors, wide range of point depths, depth discontinuities, and view-point selection for reconstruction. Experiments are provided for perspective and catadioptric cameras.

Introduction

The automatic reconstruction of photo-realistic 3D models of scenes from image sequences taken by a moving camera is still a very active field of research. Once the camera parameters of the image sequence are recovered by structure-from-motion, dense stereo and stereo merge into a single 3D model are successively applied. Currently, many 3D modeling systems exist for perspective cameras [START_REF] Nister | Automatic Dense Reconstruction from Uncalibrated Video Sequence[END_REF][START_REF] Pollefeys | Visual modeling with a handheld camera[END_REF][START_REF] Lhuillier | A quasi-dense approach to surface reconstruction from uncalibrated images[END_REF], catadioptric cameras [START_REF] Bunschoten | Robust scene reconstruction from an omnidirectional vision system[END_REF][START_REF] Lhuillier | Toward flexible 3d modeling using a catadioptric camera[END_REF] and multicamera rig [START_REF] Akbarzadeh | Towards urban 3d reconstruction from video[END_REF] (among others) sometimes with the help of additional informations such as odometry. Even if the intrinsic parameters of the camera are unknown, the involved methods are dependent on a given camera model.

Recently, it was suggested that the first sub-problem (structure-from-motion) be solved using a generic camera model and generic tools which are exploitable for any kind of camera [START_REF] Grossberg | A general imaging model and a method for finding its parameters[END_REF][START_REF] Pless | Using many cameras as one[END_REF], and the same challenge arises naturally for the complete 3D modeling process. The practical advantages would be obvious: a high ability to change one camera model for another or to mix different cameras (e.g. catadioptric camera for wide field of view and perspective camera for a few parts of the scene where higher reconstruction accuracy is needed). Many generic tools are already available for structure-from-motion: estimation of the general-ized essential matrix [START_REF] Nister | An efficient solution for the five-point relative pose problem[END_REF][START_REF] Pless | Using many cameras as one[END_REF], pose calculation [START_REF] Nister | A minimal solution to the generalized 3-point pose problem[END_REF], bundle adjustment [START_REF] Ramalingam | A generic structure-frommotion framework[END_REF][START_REF] Mouragnon | Generic and real time structure from motion[END_REF] and generic camera calibration [START_REF] Grossberg | A general imaging model and a method for finding its parameters[END_REF][START_REF] Kannala | A generic camera model and calibration method for conventional, wide-angle and fish-eye lenses[END_REF]. We have no additional contributions for this sub-problem and assume that the camera parameters are known.

Dense stereo is the second sub-problem. It is recognized that this step is very difficult in practice for uncontrolled environments. This difficulty increases in the generic context since the use of the image projection function is prohibited (this function is specific to the kind of camera). The epipolar constraint is also unavailable since the camera may be non-central. Optical flow methods [START_REF] Barron | Performance of optical flow techniques[END_REF] remain since they do not use 3D. In this second step, a hypothesis is used to obtain better results for 3D modeling: we assume that epipolar constraints are locally available in the generic images such that the standard pair-wise stereo methods [START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF] may be applied after local rectifications.

The third sub-problem is the following: once cameras and matches between image pairs are known, how can a 3D model of the scene be recovered using generic tools ? This model is a list of textured triangles in 3D which approximates the visible part of the scene where the camera has moved. We have to reconstruct 3D points, approximate them by a mesh, and deal with matching errors (false negatives and false positives), depth discontinuities, and a wide range of accuracies for reconstructed points (due to close foreground and far background, or view-point selection).

Contributions and Paper Overview

Our generic camera model is slightly different from previous ones. Previous authors [START_REF] Grossberg | A general imaging model and a method for finding its parameters[END_REF][START_REF] Ramalingam | A generic structure-frommotion framework[END_REF] model arbitrary imaging systems by a set of virtual sensing elements called raxels: a raxel is central or perspective camera with a small part of the complete view field. In our case, a raxel is reduced to a single ray (a point origin and a direction in 3D) such that the raxel center is the ray origin. We know the calibration function, which maps image pixels to rays. This function also defines the choice of all ray origins ("the ray surface choice" [START_REF] Grossberg | A general imaging model and a method for finding its parameters[END_REF]). A central camera is a special case, where all ray origins are the same point: the camera center.

Once the generic camera model is presented, Section 2 introduces a generic method to reconstruct points from 1 image matches by ray intersection, a generalization for a generic camera of virtual uncertainty [START_REF] Lhuillier | Toward flexible 3d modeling using a catadioptric camera[END_REF], and a generalization for n views of the 2-view-angle reliability [START_REF] Doubek | Reliable 3d reconstruction from a few catadioptric images[END_REF]. Both virtual uncertainty and reliability were introduced for catadioptric cameras to select reconstructed points which are retained in the final 3D model. Such selections are important if parts of the scene are to be reconstructed at very different accuracies depending on the view-point selected for reconstruction. This is also true for any other camera with a wide field of view. In this paper, both virtual uncertainty and reliability have closely related and coherent definitions, which is not the case in previous works [START_REF] Doubek | Reliable 3d reconstruction from a few catadioptric images[END_REF][START_REF] Lhuillier | Toward flexible 3d modeling using a catadioptric camera[END_REF].

Section 3 describes how to obtain a (local) 3D model for a few generic images given their corresponding camera poses, the calibration function and point correspondences. First, a reference image is chosen and segmented by a 2D mesh using gradient edges and color information. Second, points are reconstructed by our generic ray intersection. Third, 2D triangles are back-projected in 3D to fit the reconstructed points as best possible by taking into account a wide range of point depths. Virtual uncertainty is useful here to weight the minimized scores, to define the connections between triangles in 3D, and to fill holes. Finally, triangles with the worst reliability are rejected.

Section 4 provides many experiments for central cameras. Global models are obtained by combining the local models with a simple view-point selection method [START_REF] Lhuillier | Toward flexible 3d modeling using a catadioptric camera[END_REF] using our (generic) virtual uncertainty. Last, Section 5 concludes and explains what should be added for non-central cameras.

Assumptions

The proposed method involves many assumptions. First, the scene surface should be smooth enough to be approximated by a list of triangles in 3D. Second, the majority of occluding contours (and the tangent discontinuities of surfaces) should occur at gradient edges or color discontinuities in images. Third, the generic camera should not be too exotic to back-project connected image points (e.g. 2D triangles) to connected points in 3D (e.g. planar scene parts): we assume that the calibration function which maps pixels to rays in 3D is piecewise C 0 continuous with known, smooth and polygonizable discontinuities. These discontinuities occur in practice for multi-camera rig (e.g. the line between two composite images in the generic image of the stereo-rig).

Geometric Tools for a Generic Camera

This Section presents a method to reconstruct points (Section 2.1), virtual uncertainty (Section 2.2), reliability (Section 2.3) and geometric tests (Section 2.4).

The calibration function of the camera is known and maps pixels of a generic image to optical rays. An optical ray is an oriented line defined by its origin and direction. Thanks to the knowledge of camera pose in the world coordinate system, origin o and direction d (||d|| = 1) of this ray in the world coordinate system are also known and used throughout the paper.

Point Reconstruction by Ray Intersection

Once point correspondences in images, calibration and successive poses of the camera are given, 3D points of the scene should be reconstructed. The standard method to reconstruct a point is the minimization of a sum of square of reprojection errors in pixels using the Levenberg-Marquardt method [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] (LM). However, these errors cannot be used in the generic context since they require the image projection function, which is specific to the kind of camera. Many rays (o i , d i ) corresponding to observations in the i th image of the 3D point P to reconstruct should be used instead.

One solution is the reconstruction of P by minimizing the sum of squares of angles α i between vectors d i and Po i . In practice, the definition α i (P) = arccos(d i P-oi ||P-oi|| ) lead to poor LM convergence. This is not surprising: the C 2 continuity of α i is recommended for a good (quadratic and final) convergence of LM and it can be shown [1] that α i is never C 1 continuous at point P such that α i ( P) = 0.

Let R i be a rotation such that

R i d i = [0 0 1] and π the function π([x y z] ) = [x/z y/z] . Once rays (o i , d i ) are given (i ∈ {1, 2, • • • I}),
we estimate P as the minimizer of

E( P) = I i=1 ||α i ( P)|| 2 with α i ( P) = π(R i ( P -o i )). (1)
Now α i is C 2 continuous and the LM convergence is good in practice. Furthermore, ||α i (P)|| is the tangent of the angle between d i and Po i . The tangent is a good angle approximation near the expected solution where the angles are small. P is retained if it is in front of the cameras (i.e. d i (Po i ) > 0) and if E(P)/I is less than a threshold.

Virtual Covariance and Uncertainty

Assume that angle errors α i defined in Eq. 1 follow independent and identical Gaussian noise N (0 2×1 , σ 2 α I 2×2 ). Let J be the Jacobian of the function

P → [α 1 • • • α I ] .
This noise propagates to a Gaussian noise for the estimated parameter P with standard covariance matrix [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]

C(P) = σ 2 α (J(P) J(P)) -1 . (2)
An estimate of σ 2 α is obtained from residuals E(P) for all 3D reconstructed points P. Now, assume that a point P ∈ R 3 and many ray ori-

gins o i ∈ R 3 , i ∈ {1, 2, • • • I} are given. Let d i be the direction d i = P -o i ||P -o i|| .
We can solve the minimization problem defined by Eq. 1 for the new rays (o i , d i ) instead of rays (o i , d i ) and obtain the minimizer P of E with its standard covariance C(P) in Eq. 2. However,

α i (P ) = π(R i (P -o i )) = π(||P -o i ||R i d i ) = 0
and we conclude that E(P ) = 0. Since the E minimizer is unique (if P and the o i are not collinear points), we obtain P = P. Thus, the "virtual covariance matrix" of P for ray origins o i is defined by C(P ) in Eq. 2.

At this point, the expression "virtual covariance matrix" is clearer: C(P ) is the covariance matrix obtained by reconstructing P from "virtual" rays (o i , d i ) using LM, i.e. rays which are not observation rays. In the special case where P was reconstructed before by LM from other rays

(o i , d i ) corresponding to real observations in images, the corresponding standard covariance is similar to the virtual covariance if o i ≈ o i and d i ≈ d i .
Finally, the virtual uncertainty U (P ) is defined by the length of the major semi-axis of the uncertainty ellipsoid defined by C(P ) and a probability p. This ellipsoid is

∆x C -∆x ≤ X 2 3 (p), C -= 1 σ 2 α I i=1 I 3×3 -d i d i ||P -o i || 2 (3)
with X 2 3 (p) the quantile function of the X 2 distribution with 3 d.o.f, p a probability and C -the inverse [1] of C(P ).

Using notation e for the smallest eigenvalue of C -, we have

U (P ) = X 2 3 (p) e . (4)

Reliability for 3D Modeling of a Scene

A point P reconstructed from observation rays

(o i , d i ), i ∈ {1, 2,
• • • I} may be so inaccurate that the 3D model should not contain it. At first glance, we can decide that P is inaccurate for 3D modeling if U (P) is larger than a given threshold U 0 . In this Section, we introduce and justify a reliability definition R(P) which is more adequate than U (P) for thresholding.

If the generic camera is a central camera, the reconstruction is defined up to a global 3D scale and a scale change of the whole reconstruction (3D points and camera centers) implies the same scale change of the uncertainties. For a central camera, the threshold U 0 must be proportional to the scene scale to obtain a decision which is independent of the scale. This is a first reason to define

R(P) = U (P) min i ||P -o i || (5)
and decide that P is inaccurate for 3D modeling if R(P) is larger than a given threshold R 0 . We see that the permitted maximal uncertainty by condition R(P) < R 0 is proportional to the distance between point P and ray origins o i . More precisely, this inequality allows points with good accuracy (for 3D modeling of the scene) to have greater uncertainties if they are a long distance from the ray origins o i and smaller uncertainties if they are close. As a consequence, we can expect to modelize both close foreground and far background of the scene. This is the second reason for this definition of R(P). Furthermore, through this inequality it is possible to moderate the ellipsoid size (uncertainty U (P)) in comparison with the distance between ellipsoid center (the reconstructed point P) and ray origins o i . It is not difficult to prove [1] that R(P) is arbitrarily large in two cases: (1) nearly parallel d i and ( 2) large values of ||Po i ||. Case (1) occurs if all o i are collinear points and P goes toward the line of the o i . Case (2) occurs for distant point P. These cases should be avoided for 3D modeling. This is a third reason for this definition of R(P).

Geometric Tests

Let Π be the plane n X + d = 0. Once virtual covariance matrix C is defined, Mahalanobis point-to-point and point-to-plane [1] squared distances are respectively

d 2 (P 1 , P 2 ) = (P 1 -P 2 ) C -1 (P 1 )(P 1 -P 2 ) d 2 (P 1 , Π) = min P2∈Π d 2 (P 1 , P 2 ) = (n P 1 + d) 2 n C(P 1 )n . (6)
Here we introduce several tests which are systematically used by mesh operations for 3D modeling. The point-topoint neighborhood test

T (P 1 , P 2 ) is true if d 2 (P 1 , P 2 ) ≤ X 2 3 (p) and d 2 (P 2 , P 1 ) ≤ X 2 3 (p). The point-to-plane neigh- borhood test T (P 1 , Π) is true if d 2 (P 1 , Π) ≤ X 2
3 (p). The planarity test T ({P i }) is true if there is a plane Π such that all T (P i , Π) are true. In practice, Π is estimated by random samples of 3 points in the list {P i }.

These tests implicitly requires for each 3D point P i the corresponding ray origins due to the virtual covariance definition in Section 2.2. If the generic camera is central or if the points are reconstructed by LM, the ray origins are known. They are unknown in other cases, unless we apply the projection functions to P i (but this is not a generic method). More investigations are needed to estimate efficiently ray origins in the non-central case.

3D Model from Generic Images

This Section describes how to obtain a 3D model for a few generic images given their corresponding camera poses, the calibration function and image point correspondences.

First, a reference image is chosen and segmented by a 2D mesh using gradient edges and color informations (Section 3.1). Second, points are reconstructed by intersection of observation rays as described in Section 2.1. Third, 2D triangles are back-projected in 3D to fit the reconstructed points by taking into account 3D point uncertainties and depth discontinuities (Section 3.2). Last, the most unreliable parts of the resulting 2.5D mesh (Section 3.3) are rejected. Assumptions are given in Section 1.

2D Mesh

The 2D mesh in the generic reference image should satisfy many contradictory constraints: gradient edges at mesh edges, small enough mesh edges for good approximation of gradient edges, large enough mesh triangles for stable estimation of triangles in 3D and efficient rendering, uniform sampling of the field of view, and good aspect ratio for triangles. A compromise is obtained as follows.

Mesh Initialization First, a Delaunay triangulation is initialized such that the solid angles of any triangles are roughly the sames. In practice, simple checkerboards with two triangles for each rectangular cell are good enough for standard cameras like perspective, catadioptric, or stereorig. The C 0 discontinuities of the calibration function define the borders of independent 2D meshes in the reference image. Borders enforce constrained edges on the Delaunay and enforce the global shape of the checkerboards (in the catadioptric case, cell rows are concentric rings and cell columns are radial sections). The mesh resolution is defined by a mean length of cell edges equal to 8 pixels.

Gradient Edge Integration

Second, the gradient edges are integrated in the mesh by moving mesh vertices slightly and forcing mesh edges to be constrained. We have not taken into account all gradient edges since the mesh resolution has been previously fixed. So they are integrated in a best first order. A contour is a list of connected pixels which have maximum local image gradient. Its score is equal to the sum of gradient modulus for all its pixels. We pick the contour with the highest score, and find the list of closest vertices to its pixels such that the vertices have not been used before for any other contour. Then two consecutive vertices are moved slightly to approximate the contour if the part of the contour between vertex ends is a segment. Once all contours have been considered by decreasing score, a completion step is used in order to try to constrain new mesh edges if they approximate a contour in their immediate neighborhood.

Mesh Refinement Third, the 2D mesh is refined by alternating continuous improvements (move vertices to minimize a global cost combining color variance in triangles and mesh smoothness) and discrete improvements (flip edges and merge vertices to improve aspect ratio of triangles).

The continuous mesh improvement is useful for many reasons. First, few (parts of) gradient edges may be missed by the previous step, and minimizing the sum of color variances for each triangle is an other way to increase the probability that the gradient edges are on the mesh edges. Second, the gradient edge integration deformed the initial mesh only locally such that the constraint of a same solid angle for all triangles is highly violated. Minimizing the mesh smoothness (sum of squared modulus of an umbrella operator) is a way to incite incident triangles to have similar solid angles. Minimizing the mesh smoothness is also useful to improve triangle aspect ratio and regularize the minimization of color variance.

The cost function is defined by

e 2d ({p v }) = p∈t∈T ||c p - p ∈t c p |t| || 2 +λ v∈V || v ∈Nv p v -p v || 2
with T the list of mesh triangles, |t| the area of triangle t,

V the list of mesh vertices, N v the list of vertices which are connected to v by a mesh edge. Color c p at pixel p is RGB, p v is the image location of vertex v and λ is equal to 1000. The cost function is minimized using a simple descent method with vertex locations {p v } as parametrization.

All mesh vertices are allowed to move in 2D, except vertices which are incident to a constrained edge (vertices at gradient edges). The latter are only allowed to move in 1D along the detected gradient edges. This gives the priority to detected gradient edges over the minimization of color variance, which may sometimes be contradictory.

2.5D Mesh

Assume that we have I ≥ 2 camera poses and a dense list of 3D points P reconstructed by intersection of I observation rays (one for each pose) as described in Section 2.1. A 2D mesh in a reference image is also given.

First, the 2.5D mesh is initialized as a list of fully disconnected triangles in 3D. Then, this mesh is refined by alternating discrete improvements "Triangle Connection", "Hole Filling", "Triangle Removal", "Triangle Damping" and continuous improvements "Mesh Refinement". These mesh improvements are defined below thanks to the virtual covariance for the I poses (Sections 2.2 and 2.4).

At any step, the 2.5D mesh in 3D is a back-projection of the 2D mesh in the reference image. In other words, each triangle t 2d of the 2D mesh corresponds to (at most) one triangle t 3d of the 2.5D mesh with vertices v i ∈ R 3 . The t 3d vertices are parameterized by depths z i > 0 such that v i = o i + z i d i with (o i , d i ) the observation rays of t 2d vertices. Vertices in the 2D mesh may have many depths depending on current connections between triangles in 3D.

Mesh Initialization

In this step, each triangle t 2d of the 2D mesh is individually back-projected to fit the 3D points as best possible with a RANSAC procedure.

First, all 3D points reconstructed from matched pixels inside t 2d are collected in a list L t 2d . Second, planes are calculated for random samples of three 3D points taken in L t 2d . Let Π be the plane minimizing

E 2 t 2d = P∈L t 2d min{X 2 3 (p), d 2 (P, Π)} (7)
with X 2 3 (p) and d(P, Π) introduced in Eq. 3 and 6. Then, we estimate depths z i at the 3 vertices of t 2d such that o i + z i d i ∈ Π with (o i , d i ) the observation rays of these vertices. The triangle in 3D with vertices o i + z i d i is added in the 2.5D mesh if z i > 0.

Pair-Wise Triangle Connection

Triangles in 3D should be interconnected to obtain a more realistic 3D model.

Let t 3d a and t 3d b be two 3D triangles such that the associated triangles t 2d a and t 2d b in the 2D mesh have a common edge (t 2d a and t 2d b are "weakly" connected). This edge has two vertices 0 and 1 in 2D, which correspond to triangle vertices {v a 0 , v b 0 } and {v a 1 , v b 1 } in 3D. The connection between t 3d a and t 3d b is effective if the point-to-point neighborhood tests T (v a 0 , v b 0 ) and T (v a 1 , v b 1 ) defined in Section 2.4 are true.

The connection between t 3d a and t 3d b is defined as follows. Let z a i and z b i be depths such that

v a i = o i + z a i d i and v b i = o i + z b i d i with (o i , d i )
the observation rays of 2D vertices i ∈ {0, 1}. New values of z a i and z b i are set to former value of 1 2 (z a i + z b i ). Henceforth, the 2.5D mesh parameters z a i and z b i are linked by constraints z a i = z b i for further processing.

Group-Wise Triangle Connection

The "Pair-Wise Triangle Connection" above connects any triangle pair in 3D if they satisfy neighborhood conditions. Here we introduce the "Group-Wise Triangle Connection", which connects any k-group of triangles in 3D if they satisfy a planarity condition (typically k ∈ {2, 3, 4}).

A k-group of triangles in 3D is a list of k triangles t 3d j such that the corresponding triangles t 2d j are "strongly" connected in the 2D mesh. Two triangles are strongly connected if they have a common edge which is not constrained in the 2D mesh. We avoid constrained edges since they are potential surface discontinuities in 3D. Section 2.4 defines the planarity condition by T ({v i }) with {v i } the list of all triangle vertices of the k-group.

Any triangle pair {t 3d a , t 3d b } in 3D is connected as in the pair-wise case if it is included in a k-group satisfying the planarity condition and if the corresponding t 2d a , t 2d b in 2D have a common edge.

Triangle Removal A smooth surface is expected to be approximated by a list of connected triangles in 3D. If a triangle is not connected to (at least) one of its neighbors after trials of triangle connections, we have some doubt as to its quality and may decide to remove it from the 2.5D mesh. They are many reasons for fully disconnected and bad triangles in 3D: false positive matches in images (e.g. in the neighborhood of occluding contours), triangle estimations using 3D points in both close foreground and far background, too few points for reliable estimation.

Triangle Damping

The main drawback of "Triangle Removal" is the lack of triangles in scene parts which are not smooth such as tree foliages. If a triangle t 3d without connection is not removed, it may produce a major degradation of visual quality if it is very stretched in 3D in the direction d i of rays which goes across t 3d vertices. In this case, the angle θ between t 3d normal n and d i is greater than a threshold θ 0 . Thus, "Triangle Damping" reduces such degradations as follows: if θ 0 < θ, the t 3d depths z i are disturbed such that (1) the t 3d center is fixed and ( 2) n is replaced by

cos(θ 0 )d i + sin(θ 0 ) di || di|| with di = n -(n d i )d i . "
Triangle Damping" may be preferred to "Triangle Removal" to obtain more triangles in the 3D model.

Hole Filling

In our context, a hole is a connected component of triangles t 2d j in the 2D mesh without corresponding triangles t 3d j in the 2.5D mesh. "Hole Filling" is the definition of the lacking t 3d j by interpolation of depths available in the hole border. Holes are mainly due to false negative matches in low textured areas and degrade the visual quality of 3D model rendering if they are not properly filled.

The main risk is depth interpolation between foreground and background which also degrades the rendering quality, especially if foreground and background have different colors. We have the choice between strong connectivity (used in "Group-Wise Triangle Connection") and weak connectivity (used in "Pair-Wise Triangle Connection") between two triangles in the 2D mesh to define a hole as a connected component. The former is preferred to the latter which includes potential surface discontinuities at constrained edges too easily in the hole. As a consequence, the hole border is a list of edges in the 2D mesh such that (1) edges are constrained or (2) edges are not constrained and have depths at their two vertices. All 3D points corresponding to these vertices with depths are collected in a list {v i }. We also define r as the ratio between the sum of 2D lengths of edges of type (2) and the sum of 2D lengths of all border edges.

To obtain a well defined interpolation and reduce the risk of depth interpolation between foreground and background, we require that the hole border is planar thanks to the planarity condition T ({v i }) defined in Section 2.4. We also request enough 3D information at the hole border by r thresholding (0.5 < r). If T ({v i }) is true, there is a plane Π which approximates the v i and "Hole Filling" is defined as follow. Each vertex in the hole (including border) has a corresponding observation ray (o i , d i ) and a depth z i defined by o i + z i d i ∈ Π. Any hole triangle t 2d with positive z i at its vertices defines a new triangle t 3d in the 2.5D mesh.

Depth constraints are set for further processing such that these vertices have only one depth.

Mesh Refinement

The parameters of the 2.5D mesh is the list of depths z i for each triangle vertex in 3D with many constraints (equalities) between the z i . Improvements "Hole Filling" and "Pair/Group-Wise Triangle Connection" are useful to increase the rendering quality of the 3D model, but they reduce the number of independent z i and disturb the initial values of z i obtained from the 3D point cloud. The consequence is an increasing discrepancy between 3D points and the 2.5D mesh. This problem is reduced by minimizing a global cost function including a discrepancy term and a smoothness term. The smoothness term is useful to reduce noise and enforce a prior knowledge of a smooth surface on the 2.5D mesh.

The cost function to minimize is defined by

e 3d ({z i }) = t∈T E 2 t +λ {t1,t2}∈E 1 2 (|t 1 |+|t 2 |)(n t1 -n t2 ) 2
with T the list of 2D mesh triangles which has a triangle in 3D, {t, t 1 , t 2 } ⊂ T , {t 1 , t 2 } the edge between triangles t 1 and t 2 , |t| the surface (in pixels) of t, E the list of unconstrained edges in the 2D mesh, and n t the normal of the 3D triangle corresponding to the 2D triangle t. Weight λ is equal to 1 and E 2 t is defined in Eq. 7. The cost function is minimized by a descent method with depths {z i } as parametrization. Depths have a wide range due to close foreground and far background, and this should be taken into account to reduce the cost efficiently. Given a depth value z n i at iteration n of the descent method, we choose the value Step 2 merges triangles with strong conditions before step 3. Once step 3 has removed improbable and unconnected triangles, step 4 connects triangles with weaker conditions.

z n+1 i ∈ {z n i -δ i (z n i ), z n i , z n i + δ i (z n i )} which mini- mizes the partial function z i → e 3d (z i ).

Unreliable Parts

Once the 2.5D mesh is obtained, the reliability for 3D modeling (Eq. 5) allows the detection of unreliable vertices v i by thresholding such as R 0 < R(v i ). Any triangle which has an unreliable vertex is removed.

Experiments

Once camera poses and matching between image pairs are given, all experiments are done for the generic camera model restricted to central cameras (ray origins at the centre). Actually, specific pose methods [START_REF] Lhuillier | A quasi-dense approach to surface reconstruction from uncalibrated images[END_REF][START_REF] Lhuillier | Toward flexible 3d modeling using a catadioptric camera[END_REF] are prefered to generic method [START_REF] Mouragnon | Generic and real time structure from motion[END_REF] since they also estimate calibration and have more successful automatic matching. Furthermore, the dense matching method between image pairs involves local rectifications which are specific to camera model [START_REF] Lhuillier | Toward flexible 3d modeling using a catadioptric camera[END_REF].

Properties of U (P) and R(P)

In the first case (on the left of Figure 1), U (P) and R(P) are shown for two camera locations or ray origins A and B. U (P) and R(P) are defined everywhere (except on the line defined by A and B). Due to the symmetry of the problem, U and R are the same for any plane in 3D containing A and B. As expected, they increase in two cases: (1) if P goes toward line defined by A and B or (2) if P goes long away from A and B. We also see at the bottom that our reliability is very similar to the reliability given in [START_REF] Doubek | Reliable 3d reconstruction from a few catadioptric images[END_REF]: curves implicitly defined by R(P) = constant are very similar to circles defined by angle(A, P, B) = constant.

In the second case (in the middle), a camera pose C is added in the middle of A and B. The result is unexpected: there is no improvement (i.e. U or R decrease) by adding the third camera pose. In fact, the results are nearly the same. In the third case (on the right), C is moved toward the bottom. As expected, the improvement is noticeable in the neighborhood of the line defined by A and B. In these two last cases, our R definition is naturally derived from U for any numbers of views. This was not the case for the R definition of [START_REF] Doubek | Reliable 3d reconstruction from a few catadioptric images[END_REF], which was only defined for two views. radian. An image taken by this catadioptric camera is also shown.

Comparing Specific and Generic Cameras

It is wished that the virtual covariance obtained from the generic error (angle) be the same as the virtual covariance obtained from the specific error (image reprojection). We prove a simple condition [1] in image space to check this.

Assume that a point P, a specific camera model and I camera poses are given. We also consider any point X such that the image projection p i (X) is in the immediate neighborhood of p i (P) in the i-th image. The condition is

||α i (X)|| = σ α σ p ||p i (X) -p i (P)|| + o(||p i (X) -p i (P)||) (8)
with α i defined from camera center o i and direction

d i = P-o i
||P-o i || as described in Eq. 1. In other words, the projections of all circular cones (with apex at camera center) of aperture 2 radians should be circles of radius σp σα pixels. Figure 2 draws some of these projections for perspective and catadioptric cameras. In both cases, we note that the main differences between specific and generic virtual covariances occur at the borders of view fields where the circles have the largest distortions.

3D Model from Catadioptric Images

The field of view and an image taken by the catadioptric camera are also shown in Figure 2. The definition of the calibration is completed by the radii of the large and small circles: 563 and 116 pixels respectively. The image sequence has 208 images (closed turn around a church). Once the camera parameters and dense matching between image pairs are estimated, each local model is reconstructed from 3 consecutive views with σ α = 0.0011 radian and X 2 3 (0.9) = 6.25 as described in Section 3. Figure 3 shows 3D models obtained with this sequence. The local model in the first row is obtained with the reference image given in Figure 2. The most unreliable parts drawn on the left are discarded in the middle with R 0 = 0.08. A part of the ground (circular hole) and the upper part of the facade are in the blind cones defined by the small and large circles of the catadioptric images and can not be reconstructed for this reason. The second row shows views of an other local model. In both examples, we see that a successful gradient edge integration in the meshes allows sharp modeling of C 0 and C 1 depth discontinuities (in spite of the low resolution of catadioptric images).

The third row of Figure 3 shows a top view and a height map of the global model. The global model is obtained from 208 local models around the church as follows. Let U l (P) be the virtual uncertainty defined at point P with ray origins defined by the centers of cameras of a local model l. A triangle of local model l 0 is retained in the global model if U l0 (P) ≤ β min l U l (P) with β = 1.1 and P a vertex of the triangle. In other words, the triangle is retained if l 0 provides one of the best (smallest) virtual uncertainties available from all local models [START_REF] Lhuillier | Toward flexible 3d modeling using a catadioptric camera[END_REF]. This condition ignores the visibility of P in the images since U l (P) is well-defined everywhere in the generic context. In practice, the result is improved by taking into account the visibility as follows: we reset U l (P) = +∞ if P is not in the view fields of l. The global model contains 567757 triangles.

The video shows a walkthrough in the scene. The reconstruction is difficult in several parts including trees and low textured areas (e.g. street parts) which are not filled. Furthermore, the simple triangle selection above has two weakness referenced in [START_REF] Lhuillier | Toward flexible 3d modeling using a catadioptric camera[END_REF]: the model redundancy increases with depth and the self-occlusions of the surface are ignored. We have noted that the triangle selection confines the case (1) of bad reliability (Section 2.3) at the ends of image sequences, and case (1) does not occur for a closed sequence like this. Here we choose to not reject the most unreliable areas and include the far background (case (2) of bad reliability). 1 provides an estimate of the 3D noise for a few planes of the global models. The perspective camera has smaller noise than the catadioptric camera.

3D Model from Perspective Images

Conclusions

This paper presents geometric tools and results for 3D scene modeling using a generic camera model. First, virtual uncertainty and reliability are extended for generic cameras and compared with those of previous works. Second, these tools are systematically applied in the 3D model generation: fit and connect triangles in 3D, fill the holes due to matching errors, set depth resolution for 2.5D mesh optimization, re- A problem should be solved to apply the method on noncentral camera naturally: the ray origin calculation for a 3D point expressed in the camera coordinate system. Future works also include efficient matching methods in the generic context. Table 1. 3D noise for planar parts of global models. Each column provides information about a part: camera (perspective or catadioptric), number of vertices, mean distance between vertex and closest camera (cm), RMS of distances between vertex and estimated plane (cm). Vertices are selected for a part if they are projected in an ellipse of an image (white ellipses in Fig. 2 and4). Distance between two consecutive camera poses is about 30 cm.

  Virtual uncertainty is used to scale the increment δ i by δ i (z) = U (o i + zd i ) with = 0.02. Algorithm Summary Many combinations of the mesh operations above are possible and have been the subject of experiments. Our favorite strategy currently is 1. Mesh Initialization 2. apply Group-Wise Triangle Connection (k = 4), Hole Filling and Mesh Refinement alternatively 3. Triangle Removal or Triangle Damping (θ 0 = 7 20 π) 4. apply Pair-Wise Triangle Connection, Hole Filling and Mesh Refinement alternatively.

Figure 1 .

 1 Figure 1. Virtual uncertainty U and reliability R in a plane for a generic (central) camera in three cases: two camera poses on the left, three collinear poses in the middle, and three non collinear poses on the right. Camera poses are black points in this plane. Black edges are the main axes of uncertainty ellipsoids centered at some points P (their length is 2U (P). Every point P has a gray level depending on the interval where R(P) lies: [0, 1 40 [, [ 1 40 , 2 40 [, [ 2 40 , 3 40 [, [ 3 40 , 4 40 [, [ 4 40 , +∞] (darkest gray levels for largest reliabilities). On the left, we check that curves R(P) = R0 with R0 ∈ { 1 40 , 2 40 , 3 40 , 4 40 } are very similar to circles (in black). Typical values are obtained for U and R with σα = 0.001 radian and X 3 2 (0.9) = 6.25.

Figure 2 .

 2 Figure 2. Image projections of circular cones (with apex at camera center) for perspective (left) and catadioptric (middle) cameras. The former is a 35mm camera with a 30mm lens. The latter is equiangular and has a field of view of ±50 • above and below the plane orthogonal to the symmetry axis. Cone apertures are π 25

Figure 4

 4 Figure 4 shows results of the method applied to 28 (816× 1088) images taken by a perspective camera with a lateral motion. The focal length and radial distortion estimations are f = 1234 and ρ = -0.073. Local and global models are reconstructed as in the catadioptric case with σ α = 0.00059 radian, β = 1.02 and R 0 = 0.05. The global model has 129635 triangles. Table1provides an estimate of the 3D noise for a few planes of the global models. The perspective camera has smaller noise than the catadioptric camera.

Figure 3 .

 3 Figure 3. Several views of 3D models obtained with the catadioptric camera. Top and middle: views of two local models (3 consecutive poses) with rejection of unreliable triangles (R0 = 0.08), except at the top left corner. Triangle orientations are also drawn using gray levels. Bottom: top view and height map of the global model (208 poses) with rejection (R0 = 0.05).

Figure 4 .

 4 Figure 4. Left: two images (among 28) taken by the perspective camera. Right: the global 3D model.


E Mouragnon

M Lhuillier

M Dhome

F Dekeyser

P Sayd

Generic and Real-Time Structure from Motion

We introduce a generic and incremental Structure from Motion method. By generic, we mean that the proposed method is independent of any specific camera model. During the incremental 3D reconstruction, parameters of 3D points and camera poses are refined simultaneously by a generic local bundle adjustment that minimizes an angular error between rays. This method has three main advantages: it is generic, fast and accurate. The proposed method is evaluated by experiments on real data with three kinds of calibrated cameras: stereo rig, perspective and catadioptric cameras.

Introduction

The automatic estimation of scene 3D structure and camera motion from an image sequence ("Structure from Motion" or SfM) has been largely studied. Different camera models are used: pinhole, fish-eye, stereo, catadioptric, multi-cameras systems, etc. A lot of specific algorithms (i.e. specific to a given camera model) have been successfully developed and are now well known for perspective or stereo rig models [START_REF] Nister | An efficient solution to the five-point relative pose problem[END_REF][START_REF] Pollefeys | Automated reconstruction of 3D scenes from sequences of images[END_REF]. The omni-directional central (catadioptric, fish-eye) or noncentral (multi-cameras) systems that offer a larger field of view have also been widely explored [START_REF] Chang | Omni-directional structure from motion[END_REF][START_REF] Micusik | Autocalibration & 3D reconstruction with non-central catadioptric cameras[END_REF][START_REF] Pless | Using many cameras as one[END_REF]. It is a very interesting challenge to develop generic tools for SfM that are exploitable for any camera model. This way has recently been investigated with the introduction of 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] Sturm | A generic concept for camera calibration[END_REF]. In the generic camera model, pixels define image rays in camera coordinate system that can intersect or not in a unique point usually called "projection center". In recent work on generic SfM, camera motion can be estimated using generalization of the classical essential matrix [START_REF] Ramalingam | A generic structure-from-motion framework[END_REF][START_REF] Pless | Using many cameras as one[END_REF] given by 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].

A method is required for the refinement of 3D points and camera poses. The best solution for accuracy is bundle adjustment (BA) [START_REF] Triggs | Bundle adjustment -A modern synthesis[END_REF] applied on all parameters (global BA). However, it is clear that this method can not be real-time. In general, fast SfM methods or Vision-based SLAM [START_REF] Nister | Visual odometry[END_REF][START_REF] Davison | Real-time simultaneous localization and mapping with a single camera[END_REF] (Simultaneous Localization and Mapping) are less accurate than off-line methods where an optimal solution is calculated using global BA [START_REF] Triggs | Bundle adjustment -A modern synthesis[END_REF]. In this paper we present a method that makes the most of the accuracy of BA in a generic real-time application. This is possible because we developed an incremental method where not all the 3D structure is refined, but only the lastly estimated parameters at the end of the sequence. In the generic case where different cameras (pinhole, stereo, fish-eye, catadioptric...) are possible, BA is quite different from the classical one used for perspective cameras. Our generic method does not use image projections of a specific camera model but is based on back-projected rays and minimization of an angular error between rays. The first advantage is of course a high ability to change one camera model for another. The second advantage is that the method is effective if the image projection function is not explicit (as in the non-central catadioptric case) and also avoids clustering rays as in [START_REF] Ramalingam | A generic structure-from-motion framework[END_REF]. Comparison with previous works on SfM To resume, previous works are: generic but not real-time [START_REF] Ramalingam | A generic structure-from-motion framework[END_REF]. real-time but not generic [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 thanks to the use of 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 no details are given to solve this equation in common situations.

using local bundle adjustment but not generic [START_REF] Mouragnon | Real time localization and 3d reconstruction[END_REF][START_REF] Engels | Bundle adjustment rules[END_REF], and [START_REF] Engels | Bundle adjustment rules[END_REF] is not demonstrated in a real-time system with real world data as ours in this paper.

Contributions

The first contribution of our work is a generic and real-time SfM method based on an incremental 3D reconstruction and local generic bundle adjustment where an angular error is used. The second contribution is a detailed method to solve Pless Equation (in most cases, it 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]). We also compare our results with GPS ground truth and with results obtained with the most accurate (but not generic and not real-time) method available: specific global BA.

The remainder of the paper is organized as follows: Section 2 summarizes our approach and the generic camera model. The initialization method and our modified bundle adjustment are respectively explained in Section 3 and 4. Finally, experiments are presented in Section 5.

Overview of the Approach 2.1 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) where 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 unique point (camera center) whatever pixel p. In the general case, s could be any point given by calibration.

Summary

The method is based on the detection and matching of interest points (Figure 1). In each frame, Harris corners [START_REF] Harris | A combined corner and edge detector[END_REF] are detected and matched with points detected in a previous frame by computing a Zero Normalized Cross Correlation score in a near region of interest (whatever the kind of camera). The pairs with the highest scores are selected to provide a list of corresponding point pairs between the two images. To ensure a stable estimation of the 3D, a set of frames called key frames are selected. The selection criterion is based on the number of matched points between two consecutive key frames, which must be greater than a constant.

The initialization of the geometry is provided by a method based on the resolution of Pless Equation (Section 3). Then, the algorithm is incremental. For each new video frame, (1) interest points are detected and matched with those of the last key frame (2) the camera pose of the new frame is robustly estimated (Section 3.4) (3) we check if the new frame is selected as a key frame (4) if yes, new 3D points are estimated and a local bundle adjustment (Section 4) is applied. 3 Generic Initialization

The Pless Equation

Given a set of pixel correspondences between two images, the relative pose (R, t) of two cameras are 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 (q, q ′ ) such that q = d and q ′ = d ∧ s, which are convenient for this calculation. Let camera 0 be the origin of the global coordinates system and (R, t) the pose of camera 1 in this frame. The two rays must verify 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)
where [t] × is the skew symmetric cross-product matrix of the 3 × 1 vector t.

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 that we have 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 multi-camera system composed by central cameras (as the stereo rig), it means that 2D points correspondences are only made with points of the same sub-image. This is often the case in practice for two reasons: small regions of interest for reliable matching, or empty intersections between field 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 in this work that the camera motion is not a pure translation at the initialization step.

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) non-axial 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 1 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.

Camera

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

Table 1: dim(Ker(A 17 )) depends on the kind of camera.

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. It is not difficult to prove 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 3.1). 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 Ẽ. An 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. In this case, the "exact" solutions are ( R, Ẽ) ∈ V ect{(I 3×3 , 0 3×3 )} based on our assumption of "simple" matches (Section 3.1). The 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 on real data.

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 abovedescribed method and matched points are triangulated. The pose of camera 1 is estimated with 3D/2D correspondences by iterative refinement minimizing the angular error defined in Section 4.2.

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 4.2) and the overall number of counted inliers (points) includes this pose. The pose with the maximum number of inliers is then selected and another optimization is applied with all inliers.

4 Generic and Incremental Bundle Adjustment

Definitions

Bundle adjustment (BA) is the refinement of 3D points and camera poses by minimizing a cost function. The number of unknown parameters is 3 for each 3D point and 6 for each camera pose (3 for translation + for rotation). 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).

Error choice

The classical approach [START_REF] Triggs | Bundle adjustment -A modern synthesis[END_REF][START_REF] Engels | Bundle adjustment rules[END_REF] consists in the minimization of a sum of square ||ǫ i j || 2 where ǫ i j is a specific error depending on the camera model (the reprojection error in pixels). In our case, we should minimize a generic error. We define ǫ i j as the angle between the directions d i j and D i j defined above. Some experiments show that convergence of BA is bad with ǫ i j = arccos(d i j .

D i j ||D i j || ) and satisfactory with ǫ i j defined as follows (a theoretical explanation of this is given in [START_REF] Lhuillier | Automatic Scene Structure and Camera Motion using a Catadioptric System[END_REF]. We choose

ǫ 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 tangeant of the angle between d i j and D i j . The tangeant is a good approximation of the angle if it is small.

Local Generic Bundle Adjustment

In the incremental 3D reconstruction, when a new key frame I i is selected, new matched points are triangulated. Then, a stage of optimization is carried out. It 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 generic camera parameters (extrinsic parameters of key frames) and 3D points chosen for this stage i. As it is 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
Figure 2: 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.

well known that BA is very time consuming, our idea is to reduce the number of calculated parameters and avoid redundancies in computations. In our modified BA, not all the extrinsic cameras parameters are optimized but only the n last cameras parameters. Coordinates of all 3D points seen in the last n key frames are refined including new points. To bring consistency to the incremental process and ensure that new parameters are compatible with firstly estimated ones, we take account of points reprojections in the N (with N ≥ n) last frames (typically n = 3 and N = 10 are good values [START_REF] Mouragnon | Real time localization and 3d reconstruction[END_REF]). Thus, C i is the camera list {C i-n+1 . . . C i } and P i contains all the 3D points projected on cameras C i . Cost function f i is the sum of squared angular errors for all available observations in last key frames C i-N +1 . . . C i of all 3D points in P i : The incremental generic 3D reconstruction method has been tested on real data with 3 different cameras: a perspective camera, a catadioptric camera and a stereo rig. Examples of frames are available in Figure 1 and sequence characteristics in Table 2. Computation performances are reported on Table 3. In the following experiments, the trajectory obtained with our generic method is compared to GPS ground truth or global specific BA result. A rigid transformation (rotation, translation and scale factor) is applied to the trajectory as described in [START_REF] Faugeras | The representation, recognition, and locating of 3-d objects[END_REF] to fit with reference data. Then, a mean 3D error or 2D error in the horizontal plane can be measured between the generic and the reference trajectory.

f i (C i , P i ) = C k ∈{C i-N +1 ... C i } pj ∈P i ||ǫ k j || 2 . C i C i-1 C i-2 C i-3 P i C i N n

Comparison with Ground Truth (Differential GPS)

The following results are obtained with a pinhole camera embedded on an experimental vehicle equipped with a differential GPS receiver (inch precision). The vehicle trajectory is a "S" of 88 m long (Sequence 1). The calculated motion obtained with our algorithm is compared to data given by the GPS sensor and Figure 4 shows the two trajectories registration. As GPS positions are given in a metric frame we can compare camera locations and measure positioning error in meters: mean 3D error is 1.57 m and 2D error in the horizontal plane is 1.16 m. Computation time is 2 min 32 s for the whole sequence and a mean frame-rate of 6.55 f ps.

Comparison with Specific and Global Bundle Adjustment

In the two following examples, ground truth is not available. So, we compare our results with those of the best method available: a global and specific BA (all 3D parameters have been refined so as to obtain an optimal solution with a minimal reprojection error). Sequences characteristics and results are reported on Table 2. Sequence 2 is taken in an indoor environment with a hand-held pinhole camera. A very accurate result 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 3 is taken in an outdoor environment with a hand-held catadioptric camera (the 0-360 mirror with the Sony HDR-HC1E camera visible on Figure 5, DV format). The useful part of the rectified image is contained in a circle whose diameter is 458 pixels. The accuracy is also good: the mean 3D error is less than 9.2 cm for a trajectory length of (about) 40 m. The relative error is 0.23%.

Sequence 4 is taken with a stereo rig (baseline: 40 cm) in a corridor (Figure 5). The image is composed of two sub-images of 640 × 480 pix. The trajectory (20 m long) is compared to results obtained with left/right camera and global BA. The mean 3D error is 2.7/8.4 cm compared to left/right camera and the relative error is 0.13/0.42%. We have developped and experimented a generic method for the real-time Structure from Motion problem. We presented a complete process that starts with a generic initialization followed by an incremental 3D reconstruction of the scene and camera motion. The accuracy is brought by a local bundle adjustment minimizing an angular error. Experiments proved that it is easy to change one camera model for another, and promising results have been obtained on real data with three different kinds of cameras. Now, we are interested in experimenting our approach on more complex multi-camera systems.

Figure 1 :

 1 Figure 1: Feature tracks for one image of a generic camera in three cases: perspective (left), catadioptric (middle), and stereo rig (right) cameras.

Figure 3 :

 3 Figure 3: Local angular bundle adjustment when camera C i is added. Only surrounded points P i and cameras C i parameters are optimized. Nevertheless, the minimized criterion takes account of 3D points projections in the N last images.

Figure 4 :

 4 Figure 4: Left: Registration of generic vision trajectory with GPS ground truth. Continuous line represents GPS and points represent vision estimated positions. Right: 3D error (y-axis) along the trajectory (x-axis: key-frame index)

Figure 5 :

 5 Figure 5: Left: catadioptric camera and stereo rig. Middle and right: top views of 3D reconstructions for Sequence 3 (middle) and Sequence 4 (right). Trajectory in blue and 3D points in black.

Table 2 :

 2 Characteristics of video sequences.

	Sequence		Camera	#Frames #Key frames #3D Pts	#2D Pts	Traj. length
	Sequence 1		pinhole	996	66	4808	17038	88 m
	Sequence 2		pinhole	511	48	3162	11966	15 m
	Sequence 3	catadioptric	1493	132	4752	18198	40 m
	Sequence 4		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 3 :

 3 Computation times in seconds for our three cameras (detection and matching are included in Frame or Key frame times)