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

A Generic Error Model and its Application to Automatic 3D Modeling of Scenes using a Catadioptric Camera

Recently, it was suggested that structure-frommotion be solved using generic tools which are exploitable for any kind of camera. The same challenge applies for the automatic reconstruction of 3D models from image sequences, which includes structure-frommotion. This article is a new step in this direction. First, a generic error model is introduced for central cameras. Second, this error model is systematically used in the 3D modeling process. The experiments are carried out in a context which has rarely been addressed until now: the automatic 3D modeling of scenes using a catadioptric camera.

Introduction

There are two contributions in this article: the introduction of generic covariance and its application to the automatic 3D modeling of scenes using a catadioptric camera. The former and the latter are respectively summarized (and compared with previous works) in Sections 1.2 and 1.4. Section 1.1 explains why generic covariance has been introduced and Sec-tion 1.3 discusses the choice of a catadioptric camera for a visualization application. Almost all summarized results in Section 1.2 and a part of those of Section 1.4 are new materials over previous conference versions [START_REF] Lhuillier | Toward automatic 3d modeling of scenes using a generic camera model[END_REF] and [START_REF] Lhuillier | Toward flexible 3d modeling using a catadioptric camera[END_REF] of this work.

Toward Scene Modeling using Generic Tools

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] Pollefeys | Visual modeling with a hand-held camera[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 multi-camera rig [START_REF] Evers-Senne | Modelling and rendering of complex scenes with a multi-camera rig[END_REF][START_REF] Akbarzadeh | Towards urban 3d reconstruction from video[END_REF], sometimes with the help of additional information such as odometry. Even if the intrinsic parameters of the camera are unknown, the methods involved depend on a given camera model. Recently, it was suggested that the first step (SfM or Structure-from-Motion) be solved using a generic camera model and generic tools which are exploitable for any kind of camera, and the same challenge applies for the complete 3D modeling process. The practical advantage is the ability to change one camera model for another (or to mix different cameras). Many generic tools are already available for structure-from-motion: estimation of the generalized essential matrix [START_REF] Pless | Using many cameras as one[END_REF][START_REF] Li | Linear approach to motion estimation using generalized camera model[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-from-motion framework[END_REF][START_REF] Lhuillier | Effective and generic structure from motion using angular error[END_REF], 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], and even SfM system [START_REF] Mouragnon | Generic and real time structure from motion using local bundle adjustment[END_REF] (without selfcalibration).

Dense stereo is the second step. 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). Furthermore, the global epipolar constraint is unavailable since the camera may be non-central: there are no matched curves between two images such that all 3D points which project on one curve also project on the other curve. Optical flow methods [START_REF] Barron | Performance of optical flow techniques[END_REF] may be used since they do not need the global epipolar constraint. An other method is suggested in [START_REF] Lhuillier | Toward automatic 3d modeling of scenes using a generic camera model[END_REF]: assume that epipolar constraints are locally available and apply pair-wise stereo methods [START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF] after local rectifications.

The third step is the following: once cameras and matches between image pairs are known, the 3D model of the scene is reconstructed 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).

Generic covariance is introduced in [START_REF] Lhuillier | Toward automatic 3d modeling of scenes using a generic camera model[END_REF] as a tool to deal with all items of the third step.

Generic Covariance for Central Camera

Virtual covariance of a 3D point is defined for 3D modeling of scenes using a catadioptric camera [START_REF] Lhuillier | Toward flexible 3d modeling using a catadioptric camera[END_REF].

Once camera parameters are estimated by structurefrom-motion, many local 3D models along the image sequence are reconstructed by dense stereo. Then virtual covariance is used to select the points of local models which are retained in the final and global model. Generic covariance of a 3D point is defined for all central (single view point) cameras to solve the 3D modeling problem [START_REF] Lhuillier | Toward automatic 3d modeling of scenes using a generic camera model[END_REF]. Since this second covariance only depends on the 3D point and the successive positions of the camera, it is highly independent with regard to the kind of camera. Its use is also extended to other issues of 3D modeling (hole filling, surface topology decisions ...).

Note that the original name of generic covariance in [START_REF] Lhuillier | Toward automatic 3d modeling of scenes using a generic camera model[END_REF] is "virtual covariance". Here we have modified it since we think that the new name is more adequate: generic covariance does not depend on the projection function. We only assume that the camera is central.

The first theoretical contribution of the current paper is the definition of generic covariance itself. Indeed, Section 2 explains why the previous definition in [START_REF] Lhuillier | Toward automatic 3d modeling of scenes using a generic camera model[END_REF] is naive and how to obtain a thorough (mathematically sound) definition. The thorough and naive definitions provide the same final expression of generic covariance.

In short, the thorough definition for a point p and several ray origins o i is the following. First, ray directions d i corresponding to p and o i are calculated, and cost function E is defined such that E(x) is a sum of discrepancies between directions xo i and d i . Thus, p is the minimizer of E. Then, error models (random vectors) are defined for d i in the unit sphere. There is a first-order error propagation of these errors to the minimizer of E, and the generic covariance of p is defined as the covariance of this minimizer.

One of the interests of generic covariance is its simplicity. It only depends on p, o i and a scale factor σ α .

The second theoretical contribution of the paper is a list of properties of generic covariance. Section 3 provides the links (1) between unit sphere error and image error, (2) between virtual covariance [START_REF] Lhuillier | Toward flexible 3d modeling using a catadioptric camera[END_REF] and generic covariance, (3) between σ α and ray intersection problems from image points. In case [START_REF] Akbarzadeh | Towards urban 3d reconstruction from video[END_REF], we give a simple criterion to check if both covariances are equal. In case [START_REF] Bunschoten | Robust scene reconstruction from an omnidirectional vision system[END_REF], we explain how to estimate σ α . Section 4 provides asymptotic properties of the uncertainty derived from generic covariance. Here we call "uncertainty" the length of the major semiaxis of the ellipsoid defined by generic covariance and a given probability. We show that uncertainty of a point p increases as the square of the distance between point p and the ray origins {o i }. This is a generalization for all central cameras of a well known property for rectified stereo-rig: depth z = Bf d has uncertainty σ z = z 2 Bf σ d using B=baseline, f =focal length, d=image disparity and σ d =disparity uncertainty.

There is an other asymptotic property observed in experiments [START_REF] Lhuillier | Toward automatic 3d modeling of scenes using a generic camera model[END_REF] and proven in Section 4. In the two view case, the ratio of generic uncertainty by distance between p and {o 1 , o 2 } is the inverse of apical angle (up to a constant scale). The apical angle defined by a point p and two ray origins o 1 and o 2 is the angle between po 1 and po 2 . There is also a generalization in the multi-view case. The ratio between our uncertainty and distance is called "reliability".

This relation between apical angle and generic uncertainty has consequences on many previous works. On the one hand, reliability (derived from generic covariance) is thresholded in [START_REF] Lhuillier | Toward automatic 3d modeling of scenes using a generic camera model[END_REF] to reject reconstructed points in 3D models which are considered "unreliable". On the other hand, apical angles are also thresholded to reject points [START_REF] Doubek | Reliable 3d reconstruction from a few catadioptric images[END_REF] for the same reason. A criterion for video sub-sampling (used to stabilize structure-from-motion) is also based on apical angles [START_REF] Torii | Measuring camera translation by the dominant apical angle[END_REF]. Now we see the link between these methods: heuristic methods based on apical angles are roughly equivalent to generic covariance based methods. The former are two-view based and need recipes for more than two views. The latter are intrinsically multiview.

Catadioptric Camera for Visualization Application

Our target application is the automatic reconstruction of photo-realistic 3D models for walkthroughs in a complex scene. This is a long-term research problem in Computer Vision and Graphics. A minimal requirement for interactive walkthrough is the scene rendering in any view direction around the horizon-tal plane, when the viewer moves along the ground. This suggests a wide field of view for the given images, for which many kinds of camera are possible [START_REF] Daniilidis | The page of omnidirectional vision[END_REF]: catadioptric cameras, fish-eyes, or systems of multicameras pointing in many directions. Since we would like to capture any scene where a pedestrian can go, the hardware involved should be hand-held/headheld and not cumbersome. A catadioptric camera is a good candidate for all these constraints. The main drawback of this choice is the low resolution compared with a perspective camera for a given field of view. We compensate for this problem using still image sequences taken with an equiangular catadioptric camera. Still images are preferred to video images because of their better quality (resolution, noise). An equiangular camera has also been selected from among other catadioptric cameras, since it is designed to spread the resolution of view field well throughout the whole image. Today, such a catadioptric camera can be purchased on the web for a few hundred Euros using adequate mirrors [1].

These two choices (still images and equiangular camera) mainly have two consequences. First, a still image sequence requires some effort and patience: the user should alternate a step in the scene and a (non blurred) shot by pressing a button. Second, an equiangular catadioptric camera is not a central camera. A non-central model complicates involved methods, since the back-projected rays do not intersect a single point in space as the perspective model.

This paper shows that the results obtained are worthwhile with such a setup, and that a central approximation of equiangular camera is sufficient in many cases.

Scene Modeling using Catadioptric Camera

Our 3D modeling methods are essentially generic for central cameras thanks to the systematic use of generic covariance. We integrate them in a fully automatic reconstruction system for catadioptric image sequences and experiment on hundreds of images without precise calibration knowledge. The overall system and experiments themselves are practical contributions of this work. Indeed, fully automatic 3D modeling from catadioptric images has rarely been addressed until now (although this is a long-standing problem for perspective images). All previous catadioptric approaches have been limited to dense reconstruction for a few view points or use accurate calibration of the camera. Experiments in Section 7 include outdoor scene reconstructions (only indoor examples are provided in previous works [START_REF] Kang | 3-d scene data recovery using omnidirectional multibaseline stereo[END_REF][START_REF] Bunschoten | Robust scene reconstruction from an omnidirectional vision system[END_REF][START_REF] Fleck | Omnidirectional 3d modeling on a mobile robot using graph cuts[END_REF]). Now, the system is summarized. Structure-from-motion (SfM) is the first step. Although the principles are well known, optimal (including global bundle adjustment) and robust SfM systems are not so common if the only given data is a long image sequence acquired by a general catadioptric camera and imprecise knowledge about the calibration. Previous catadioptric SfM systems are [START_REF] Strelow | Extending shape-from-motion estima-tion to noncentral omnidirectional camera[END_REF][START_REF] Micusik | Structure from motion with wide circular field of view cameras[END_REF][START_REF] Lhuillier | Automatic scene structure and camera motion using a catadioptric system[END_REF]. The method published in [START_REF] Lhuillier | Automatic scene structure and camera motion using a catadioptric system[END_REF] is used in this work. This method is preferred to the generic method [START_REF] Mouragnon | Generic and real time structure from motion using local bundle adjustment[END_REF] since it also estimates camera calibration and it has more successful automatic matching.

Here we use the central camera model with a general radial function, which is estimated from an approximate knowledge of the two view field angles provided by the mirror manufacturer. Details are omitted here since they have been published before.

Dense stereo is the second step. Catadioptric images are reprojected on virtual cubes, then dense stereo is applied on parallel faces such that conjugated epipolar lines are parallel. This is a particular case of the dense stereo generic scheme referenced in Section 1.1 and [START_REF] Lhuillier | Toward automatic 3d modeling of scenes using a generic camera model[END_REF]. There are at least two reasons for using virtual cubes [START_REF] Lhuillier | Toward flexible 3d modeling using a catadioptric camera[END_REF] instead of virtual cylinders [START_REF] Kang | 3-d scene data recovery using omnidirectional multibaseline stereo[END_REF][START_REF] Bunschoten | Robust scene reconstruction from an omnidirectional vision system[END_REF][START_REF] Fleck | Omnidirectional 3d modeling on a mobile robot using graph cuts[END_REF] or even the original catadioptric images. First, a large choice of dense stereo methods is available [START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF]. Second, this facilitates the 3D reconstruction of the scene ground by pre-rectification. Here we do not focus on a specific stereo method. For convenience, we use the quasi-dense propagation method [START_REF] Lhuillier | Match propagation for image-based modeling and rendering[END_REF], but other methods are possible: [START_REF] Kang | 3-d scene data recovery using omnidirectional multibaseline stereo[END_REF][START_REF] Bunschoten | Robust scene reconstruction from an omnidirectional vision system[END_REF] use multi-baseline stereo inspired by [START_REF] Okutomi | A multiple-baseline stereo[END_REF], [START_REF] Fleck | Omnidirectional 3d modeling on a mobile robot using graph cuts[END_REF] uses graph-cut method [START_REF] Kolmogorov | Computing visual correspondances with occlusions using graphcuts[END_REF] followed by postprocessings.

The third step is the 3D model generation from camera poses and image matches using generic covariance. Section 5 describes how to obtain a local (view-centered) 3D model for a few images and discusses limitations. First, a reference image is seg-mented by a 2D mesh using gradient edges and color information. Second, points are reconstructed from image matches by ray intersection. Third, 2D triangles are back-projected in 3D to fit the reconstructed points. Generic covariance is useful here to weight the minimized scores, to define the connections between triangles in 3D, to fill holes, and to reject unreliable triangles. Once many local 3D models have been reconstructed along the image sequence, generic covariance is once again used to obtain the final and global 3D model by view point selection and redundancy reduction (Section 6). Note that view point selection is a key issue [START_REF] Lhuillier | Toward flexible 3d modeling using a catadioptric camera[END_REF]: a 3D point of the scene may be reconstructed in several local models at very different accuracies, and local models with the worst accuracies must not be used to reconstruct this point.

Definition of Generic Covariance

This Section provides two generic covariance definitions for central camera. The first one [START_REF] Lhuillier | Toward automatic 3d modeling of scenes using a generic camera model[END_REF] is given in Section 2.2. Then Section 2.3 explains (1) why this definition is "naive" and (2) how to obtain a "thorough" definition. Last, the three steps of the thorough definition are set out in Sections 2.4, 2.5 and 2.6.

Notations

Several notations are used throughout the paper. Different fonts are used for reals (eg. 0, 1, a, σ), vectors (eg. b), matrices (eg. C) and functions (eg. F (x), g(x)). The real vector space of dimension k is R k . The identity matrix of dimension k is I k . Let k be the vector 0 0 1 T and 0 k be the null vector

of R k . The angle between two vectors a, b ∈ R 3 is (a, b) ∈ [0, π] and the cross product is a ∧ b. The unit sphere (not ball) of R 3 is S 2 . The set of rotations of R 3 is SO(3). Let π be the function π( x y z T ) = x z y z T . (1)
Notation z ∼ N (z, C z ) means that z is a Gaussian vector which has mean z and covariance C z .

Naive Definition

Let o i ∈ R 3 , i ∈ {1, 2, • • • , I} be many ray origins and a point p ∈ R 3 \ {o i } such that p and {o i } are not collinear. We introduce directions d i ∈ S 2 and choose rotations R i ∈ SO(3) such that

d i = p -o i ||p -o i || and R i d i = k. ( 2
)
Point p is a minimizer of the cost function

E(x) = I i=1 ||α i (x)|| 2 with α i (x) = π(R i (x -o i )) (3) since E(p) = 0. Using notation z = R i (x -o i ) = x y z T , we have ||π(z)|| 2 = x 2 + y 2 z 2 = tan 2 (k, z) = tan 2 (d i , x -o i ). ( 4
)
Thus, E(x) is the sum of squares of tangents of angles (d i , xo i ). Note that p is the unique minimizer of E since p and {o i } are not collinear. Let σ α > 0. We assume [START_REF] Lhuillier | Toward automatic 3d modeling of scenes using a generic camera model[END_REF] that the angle errors α i follow independent, isotropic and identical Gaussian errors N (0 2 , σ 2 α I 2 ). Let J be the Jacobian of the function

x → α T 1 (x) • • • α T I (x) T .
There is a first-order error propagation from the α i to the minimizer p of E: p follows Gaussian error with covariance matrix

C(p) = σ 2 α (J(p) T J(p)) -1 . (5)
The next step is to calculate C(p). Let J π and J αi be the Jacobians of π and α i . The Chain rule provides

J αi (p) = J π (||p -o i ||k)R i = 1 ||p -o i || AR i (6)
using A = 1 0 0 0 1 0 . Thanks to Eqs. 5 and 6,

C -1 (p) = 1 σ 2 α J(p) T J(p) = 1 σ 2 α I i=1 1 ||p -o i || 2 R T i A T AR i = 1 σ 2 α I i=1 1 ||p -o i || 2 R T i (I 3 -kk T )R i . (7) Since R T i R i = I 3 and R T i k = d i ,
we obtain a simple expression of the generic covariance matrix:

C(p) = σ 2 α ( I i=1 I 3 -d i d ⊤ i ||p -o i || 2 ) -1 . (8)
Note that E and C(p) do not depend on the R i choice, but α i and J do.

Against the Naive Definition

The first-order error propagation (for Gaussian vector) is formulated like this [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. Assume that φ is a C 1 continuous function, J φ is the Jacobian of φ, and y ∼ N (y 0 , C y ). Up to the first order, y propagates to

φ(y) ∼ N (φ(y 0 ), C φ ) with C φ = J φ (y 0 )C y J T φ (y 0 ). (9)
In the naive definition, Eq. 5 is obtained by propagation as if there were a function φ which maps value of α T 1

• • • α T I
T to the minimizer of E. However, the value of E is fixed by the values of all α i . In this case, the minimizer of E and φ itself are not well defined.

A correct use of propagation is the following:

1. Choose an error model for vectors {d i , o i } 2. Estimate Jacobian of function φ which maps model y of (d

1 , o 1 , • • • , d I , o I ) to the minimizer p of E 3.
Define the generic covariance C(p) as the covariance of φ(y) using Eq. 9.

Henceforth, the subject of Section 2 is the development of these steps. Section 2.4 presents an error model choice (step 1) such that Eq. 5 is still correct. So Eq. 8 is too. The main part is the estimation of J φ (step 2) in Section 2.5. Section 2.6 describes the so-called "first-order error propagation" (step 3).

Error Model of Ray Directions d i and Origins o i

An ideal error model d * i of d i modelizes perturbations of d i on the unit sphere S 2 since d i ∈ S 2 . In this article, we simplify the problem by using a Gaussian error model d * i which approximately modelizes perturbations on S 2 : d * i modelizes isotropic perturbations on T , the tangent plane on S 2 at d i .

Let {a, b} be an orthonormal basis of T , σ α > 0,

n i ∼ N (0 2 , I 2 ) and d * i = d i + σ α a b n i . Then, propagation from n i to d * i implies d * i ∼ N (d i , C i ) and C i = σ 2 α a b I 2 a b T = σ 2 α a b d i (I 3 -kk T ) a b d i T = σ 2 α (I 3 -d i d T i ). ( 10
)
Last we assume that the d * i are independent with the same scale σ α and the o i have no errors.

Estimation of Jacobian

J φ Let R i be a C 2 continuous function from a neighbor- hood of d i into SO(3) such that R i (x)x = ||x||k. (11)
Now we define the cost function

E * (x) = I i=1 ||α * i (x)|| 2 , α * i (x) = π(R i (d * i )(x -o i )) (12)
and its minimizer p * . The goal of Section 2.5 is the estimation of Jacobian J φ of function φ, which maps

(d * 1 , d * 2 , • • • , d * I ) to p * .
Note that E * is the sum of squares of the tangent of (d * i , x-o i ). Thus E * and p * do not depend on the R i choice, and we have E * = E and p * = p if ∀i, d * i = d i . Furthermore, φ does not depend on error model of ray origins since this model is not defined in our case (Section 2.4). Lemma 1 provides the existence of R i .

Lemma 1 Assume that a, b ∈ S 2 . There is a C 2 continuous function x → R x from the half space {x ∈ R 3 , a T x > 0} into SO(3) such that R x x = ||x||b.
Proof First, assume a T b = 0 and a T x > 0. Thus x is not parallel to b and x ∧ b = 0. Let n = x∧b ||x∧b|| . Since x and ||x||b have the same norm and are both orthogonal to n, the rotation R x defined by axis direction n and angle (x, b) is such that R

x x = ||x||b. Function x → R x is C 2 continuous since functions x → (n, (x, b)) and (n, (x, b)) → R x are C 2 continu- ous.
Second, assume a T b = 0. Let R 0 be a rotation such that a T R 0 b = 0. We use the previous scheme to obtain a function R

′ x such that R ′ x x = ||x||R 0 b with x in the half space a T x > 0. So function R x : x → R T 0 R ′ x
is such that R x x = ||x||b for all x in the half space a T x > 0. Thanks to this Lemma, we know that

E * is well defined if d * i is in the half space H i = {x ∈ R 3 , d T i x > 0}.
The error model of Section 2.4 meets d * i ∈ H i . Now, we rewrite Eq. 12 in a different form to simplify the J φ estimation. Let

d * i ∈ H i , x ∈ R 3 \ {o i }, y T = (d * 1 ) T • • • (d * I ) T (13)
and F be the function

F (x, y) = α * 1 (x) T • • • α * I (x) T T . (14)
We obtain

E * (x) = ||F (x, y)|| 2 and φ(y) = arg min x E * (x). ( 15
)
Note that the estimation of J φ is not the same as the standard estimation encountered in 3D Vision. E.g. for bundle adjustment, we should estimate J φ such that

φ(y) = arg min x ||F (x) -y|| 2 (16)
with x the 3D parameters (cameras poses and scene points), y the points detected in images, and F the projection functions. Then we have

J φ ≈ ( ∂F ∂x T ∂F ∂x ) -1 ∂F ∂x T . ( 17
)
Our case is more general and is solved by Lemma 2.

Lemma 2 Let F be a C 2 continuous function with full rank Jacobian ∂F ∂x . There is a C 1 continuous function

φ(y) = arg min x ||F (x, y)|| 2 (18)
such that

J φ ≈ -( ∂F ∂x T ∂F ∂x ) -1 ∂F ∂x T ∂F ∂y ( 19
)
with derivatives of F taken at (φ(y), y). This approximation is exact equality if F (φ(y), y) = 0.

Proof Lemma 2 is a particular case of Proposition 6.1 in [START_REF] Faugeras | Geometry of Multiple Images[END_REF], which also asserts that function φ locally exists and is C 1 continuous. Notations F k , x i and y j are the k-th, i-th and j-th coordinates of function F and vectors x and y. The coefficient at i-th row and j-th column of matrix M is M i,j . Second-order partial derivatives of

f (x, y) = 1 2 ||F (x, y)|| 2 (20)
are

∂ 2 f ∂x i ∂y j = k { ∂F k ∂x i ∂F k ∂y j + ∂ 2 F k ∂x i ∂y j F k }. (21)
The Gauss-Newton approximation of this equation is

∂ 2 f ∂x i ∂y j ≈ k ∂F k ∂x i ∂F k ∂y j = ( ∂F ∂x T ∂F ∂y ) i,j . (22)
A similar approximation is

∂ 2 f ∂x i ∂x i ′ ≈ k ∂F k ∂x i ∂F k ∂x i ′ = ( ∂F ∂x T ∂F ∂x ) i,i ′ . ( 23
)
These approximations are exact equalities if

F (x, y) = 0. Since φ(y) is minimizer of x → f (x, y), we have ∀i, ∂f ∂x i (φ(y), y) = 0. ( 24
)
Using Eqs. 24, 22 and 23, we deduce that ∀i, ∀j,

0 = ∂ ∂y j (y → ∂f ∂x i (φ(y), y)) = ∂ 2 f ∂x i ∂y j + i ′ ( ∂ 2 f ∂x i ∂x i ′ ) ∂φ i ′ ∂y j ≈ ( ∂F ∂x T ∂F ∂y ) i,j + i ′ ( ∂F ∂x T ∂F ∂x ) i,i ′ ( ∂φ ∂y ) i ′ ,j ≈ ( ∂F ∂x T ∂F ∂y + ( ∂F ∂x T ∂F ∂x ) ∂φ ∂y ) i,j . (25)
Thanks to the full rank of ∂F ∂x , matrix ( ∂F ∂x T ∂F ∂x ) is invertible and we obtain the result.

Last, the following Lemma explicits ∂F ∂x and ∂F ∂y .

Lemma 3 Let y T 0 = d 1 T • • • d I T , R i = R i (d i )
and assume d i = p-oi ||p-oi|| . We have

∂F ∂x =    A 1 . . . A I    and ∂F ∂y =    B 1 0 0 0 . . . 0 0 0 B I    (26)
where

A i = AR i ||p -o i || , B i = -AR i , A = 1 0 0 0 1 0 . (27)
The derivatives of F are taken at (x, y) = (p, y 0 ).

Proof The block-wise structures of F derivatives result from the definition of F . First we calculate

∂α * i ∂x at (x, d * i ) = (p, d i )
. Thanks to Eq. 11 and

d i = p-oi ||p-oi|| , R i (d i )(p -o i ) = ||p -o i ||k. ( 28
)
We also need the Jacobian of π (π is defined in Eq. 1)

J π ( x y z T ) = 1 z 0 -x z 2 0 1 z -y z 2 . ( 29
)
Then we apply the Chain rule (as in Eq. 6)

A i = ∂α * i ∂x (p, d i ) = J π (||p -o i ||k)R i = AR i ||p -o i || . (30)
Second we calculate

∂α * i ∂d * i at (x, d * i ) = (p, d i ). Using d i = p-oi ||p-oi|| , we have α * i (p) = π(R i (d * i )(p -o i )) = π(R i (d * i )d i ). ( 31
)
The Chain rule and

R i (d i )d i = k provide ∂α * i ∂d * i (p, d i ) = J π (k) ∂(R i (d * i )d i ) ∂d * i (d i ) = A ∂(R i (x)d i ) ∂x (d i ). (32)
Fortunately, an explicit expression of R i is not needed to calculate

∂α * i ∂d * i (p, d i ):
we derivate Eq. 11 with respect to parameter a ∈ {x, y, z} of x = x y z

T ∂(R i (x)) ∂a x + R i (x) ∂x ∂a = k ∂||x|| ∂a (33)
and we obtain at point

x = d i ∂(R i (x)d i ) ∂x (d i ) + R i = kd T i . (34)
We deduce from Eqs. 32 and 34 that

B i = ∂α * i ∂d * i (p, d i ) = A(kd T i -R i ) = -AR i (35)
Thanks to Lemma 3, ∂F ∂x has full rank if p and {o i } are not collinear. Then Lemma 2 is used: φ locally exists, it is C 1 continuous, and its Jacobian at (p, y 0 ) is easy to estimate thanks to Eqs. 26 and 27.

First-Order Error Propagation

The last step of the generic covariance definition is the use of first-order error propagation (Eq. 9) with

y T = (d * 1 ) T • • • (d * I )
T , φ of Section 2.5 and its Jacobian (Eq. 19):

C φ = ( ∂F ∂x T ∂F ∂x ) -1 ∂F ∂x T ∂F ∂y C y ∂F ∂y T ∂F ∂x ( ∂F ∂x T ∂F ∂x ) -1 . (36) at (x, y) = (φ(y 0 ), y 0 ) = (p, d 1 T • • • d I T T ). Thanks to the d * i definition in Section 2.4, C y is a block-wise diagonal matrix with diagonal blocks equal to C i . Then Lemma 3 implies that ∂F ∂y C y ∂F ∂y
T is a block-diagonal matrix and its diagonal blocks are

B i C i B T i = σ 2 α AR i (I 3 -d i d T i )R T i A T = σ 2 α A(I 3 -kk T )A T = σ 2 α I 2 (37) if d i = p-oi ||p-oi|| . Thus, the generic covariance C(p) is C φ = σ 2 α ( ∂F ∂x T ∂F ∂x ) -1 = σ 2 α (J(p) T J(p)) -1 . (38)
This result is the same as that of the naive definition.

Observing Eq. 37, we note that alternative errors d * i are possible to obtain the same C(p): conditions are

C i = σ 2 α I 3 + λ i d i d T i , λ i ∈ R (d i is symmetry axis of C i ). The alternative errors include non Gaussian errors d * i ∈ S 2
, without the tangent plane approximation of Section 2.4.

Properties of Generic Covariance

First, Section 3.1 investigates how the S 2 error (introduced in Section 2.4) is propagated to image space. Then Section 3.2 provides a criterion to check if the generic covariance propagated from S 2 and the covariance propagated from image space are the same.

In both cases, the input error is isotropic (in S 2 or image) and has uniform scale (σ α or σ p ). Last, Section 3.3 explains how to estimate σ α .

Image Error and S 2 Error

The error on the unit sphere S 2 propagates to image error by the projection function of the central camera. This propagation is described in Lemma 4.

Lemma 4 Let p be the projection function of the camera, and d ∈ S 2 . Assume that p is C 1 continuous with full rank Jacobian J p . Up to the first order, the S 2 error

d * ∼ N (d, C) with C = σ 2 α (I 3 -dd T ). ( 39
)
propagates to image error 40) Furthermore, the image projection of the circular cone with infinitesimal aperture 2ǫ radians, apex 0, and axis d is in the ellipse

p(d * ) ∼ N (p(d), C p ) with C p = σ 2 α J p (d)J T p (d). (
{m ∈ R 2 , (m -p(d)) T ( ǫ 2 σ 2 α C p ) -1 (m -p(d)) = 1}. ( 41
)
Proof Point zd defined by direction d and z > 0 is such that p(zd)p(d) = 0 since 0 is the camera center. Using z and d derivatives of this equation, we obtain

J p (zd)d = 0 and zJ p (zd) -J p (d) = 0. ( 42
)
Thus Eq. 39 propagates to p(d * ) ∼ N (p(d), C p ) with

C p = σ 2 α J p (d)(I 3 -dd T )J T p (d) = σ 2 α J p (d)J T p (d). ( 43
)
Now, we use SVD and introduce several notations:

J p (d) = UDV T , R = V d , x = R x y z T . (44) Matrix R is a rotation since d ∈ S 2 and J p (d)d = 0. Thus, ǫ 2 z 2 = x 2 + y 2 iff
x is in the circular cone with apex 0, axis d and aperture 2 arctan(ǫ) ≈ 2ǫ radians.

The linear Taylor expansion of p at d and Eq. 44 imply

p(x) -p(d) = p( 1 z x) -p(d) ≈ J p (d)( 1 z x -d) ≈ J p (d)R x z y z 0 T = 1 z UD x y . (45)
Note that Eq. 45 approximation is correct for x in the cone with small enough ǫ since

|| 1 z x -d|| 2 = x 2 +y 2 z 2 = ǫ 2 .
Furthermore, D > 0 since J p has full rank. Now, the SVD of J p (d), invertible D, Eqs. 43 and 45 provide

(p(x) -p(d)) T ( ǫ 2 σ 2 α C p ) -1 (p(x) -p(d)) ≈ 1 z 2 x y DU T (ǫ 2 UDV T VDU T ) -1 UD x y T ≈ x 2 + y 2 z 2 ǫ 2 . ( 46
)
We obtain the last result since ǫ 2 z 2 = x 2 + y 2 iff x is in the circular cone.

Lemma 4 also provides a method to visualize the covariance matrix C p of the propagated image error p(d * ): uncertainty ellipse of image error p(d * ) is distortion ellipse of projection p centered at p(d). We must remember that the distortion ellipse (or Tissot's indicatrix [START_REF] Tissot | Tissot's indicatrix[END_REF]) of p is the image projection by p of circular cone with infinitesimal aperture 2ǫ radians, apex 0 and axis d.

Thus, there is a convenient and visual test to check if image error p(d * ) is "standard" (isotropic and uniform in the whole image): the more ellipses are circles with the same radius, the more standard the image error. Note that unavoidable distortions occur by local map from sphere into plane.

Generic and Virtual Covariances in 3D

On the one hand, the generic covariance C(p) is defined by error propagation from the unit sphere S 2 (Section 2). On the other hand, the virtual covariance C v (p) is defined by error propagation from the image [START_REF] Lhuillier | Toward flexible 3d modeling using a catadioptric camera[END_REF]. Both are covariances for a 3D point p.

The former (the latter, respectively) assumes that the error amplitude in S 2 (in the image, respectively) is uniform. The latter is specific to the projection function p of the camera. The following Lemma provides the condition to obtain

C v = C. Lemma 5 Let p be the projection function of the camera, R 0 i ∈ SO(3), o i ∈ R 3 and p i (p) = p(R 0 i (p -o i )). ( 47
)
Assume that p is C 1 continuous. Let J pi be the Jacobian of p → p i (p), σ p > 0 and

C v (p) = σ 2 p ( I i=1 J T pi (p)J pi (p)) -1 . ( 48
)
The S 2 error is

∀d ∈ S 2 , d * ∼ N (d, C) with C = σ 2 α (I 3 -dd T ). (49) If all p(d * ) have the same covariance σ 2 p I 2 , C v = C. Proof According to Lemma 4, p(d * ) ∼ N (p(d), C p ) with C p = σ 2 α J p (d)J T p (d). Furthermore, the SVD J p (d) = UDV T and C p = σ 2 p I 2 imply σ 2 p I 2 = C p = σ 2 α UDV T VDU T ⇒ D = σ p σ α I 2 . (50)
Thus, there are U, V and R = V d such that

J p (d) = σ p σ α UV T and U T U = I 2 , V T d = 0 and VV T = R(I 3 -kk T )R T = I 3 -dd T . (51) Let d i be such that p-o i = ||p-o i ||d i .
Using Chain Rule, Eqs 42 and 51, the Jacobian of p i at p is

J pi (p) = J p (R 0 i (p -o i ))R 0 i = 1 ||p -o i || J p (R 0 i d i )R 0 i = 1 ||p -o i || σ p σ α U i V T i R 0 i ( 52
)
with

U T i U i = I 2 , V i V T i = I 3 -R 0 i d i (R 0 i d i ) T . (53)
Now, Eqs. 48, 52 and 53 imply

C -1 v (p) = 1 σ 2 p I i=1 J T pi (p)J pi (p) = 1 σ 2 p I i=1 1 ||p -o i || 2 σ 2 p σ 2 α R 0 i T V i U T i U i V T i R 0 i = 1 σ 2 α I i=1 1 ||p -o i || 2 (I 3 -d i d T i ) = C -1 (p). ( 54
)
Lemma 5 asserts that generic and virtual covariances of a 3D point are the same if the image error is standard, i.e. if all p(d * ) have the same isotropic covariance. Thus, the visual test provided at the end of Section 3.1 can be used to check if generic and virtual covariances are the same.

In practice, perspective cameras with moderated field of view have distortion ellipses which are similar to circles with the same radius. According to Section 3.1, these cameras propagate the S 2 error to the standard image error. According to Lemma 5, virtual and generic covariances are similar for perspective cameras.

Ray Intersection from Points in Images

Assume that we have a ray intersection problem: there are noisy and matched points in I images such that camera calibration and camera poses are known, and the corresponding 3D point p should be reconstructed. This problem may be solved as follows:

1. Calculate the ray directions d i and origins o i corresponding to image points (in world coordinates).

2. Choose rotations R i such that R i d i = k and define the angle cost function E(x) as in Eq. 3.

Estimate p as the minimizer of E(x).

Here we can not define a covariance for p using Eq. 38 since it requires d i = p-oi ||p-oi|| , which is wrong due to reconstruction in the presence of image noise. However, Eq. 36 does not need d i = p-oi ||p-oi|| and we use it to define covariance C r (p). Now a 3D point p is reconstructed, we can reset d i using d i = p-oi ||p-oi|| and define covariance C(p) using Eq. 8. Both covariances C r and C are obtained, the former comes from points detected in images and the latter comes from a reconstructed point in space.

Assume that the image noise is low, the mapping from image point to ray direction is continuous, and the mapping from ray directions to covariance is continuous by Eq. 36. Then C r is approximated by C. This approximation is done in all our experiments.

Last we present a Lemma which is used in experiments to estimate σ α from several ray intersection problems as defined above.

Lemma 6

The mean (expected value) of random variable E * (φ(y)) is

E(E * (φ(y))) ≈ (2I -3)σ 2 α . (55)
Proof Here we need new notations D(y) = F (φ(y), y) and

P = ∂F ∂x ( ∂F ∂x T ∂F ∂x ) -1 ∂F ∂x T
. Note that P 2 = P = P T . Furthermore, Section 2.5 and Eq. 37 provide several relations which are useful for the current proof:

y T 0 = d T 1 • • • d T I φ(y 0 ) = p, D(y 0 ) = F (p, y 0 ) = 0 ∂F ∂x J φ = -P ∂F ∂y , ∂F ∂y C y ∂F ∂y T = σ 2 α I 2I . (56)
Partial derivatives of F are taken at (x, y) = (p, y 0 ). Now, a linear Taylor expansion of D is

D(y) ≈ D(y 0 ) + ∂D ∂y (y 0 ).(y -y 0 ) ≈ ( ∂F ∂x (p, y 0 )J φ (y 0 ) + ∂F ∂y (p, y 0 )).(y -y 0 ) ≈ (I 2I -P) ∂F ∂y (p, y 0 ).(y -y 0 ). ( 57
)
Using first-order error propagation of y ∼ N (y 0 , C y ), we have D(y) ∼ N (0 2I , C D ) such that

C D = (I 2I -P) ∂F ∂y C y ∂F ∂y T (I 2I -P) T = σ 2 α (I 2I -P -P T + PP T ) = σ 2 α (I 2I -P). ( 58
)
Now the expected value of E * (φ(y)) is

E(E * (φ(y))) = E(D(y) T D(y)) = E(tr(D(y)D(y) T )) = tr(E(D(y)D(y) T )) ≈ tr(C D ) = 2Iσ 2 α -σ 2 α tr(P) = 2Iσ 2 α -σ 2 α tr(( ∂F ∂x T ∂F ∂x ) -1 ∂F ∂x T ∂F ∂x ) = σ 2 α (2I -3). ( 59
)
This lemma is used as follows. We reconstruct J points using the method described at the beginning of Section 3.3 for a same number of I images. Let E j be the final value of the minimized score of the j-th point. Thanks to Eq. 55, we estimate σ α using

(2I -3)σ 2 α ≈ E(E * (φ(y))) ≈ 1 J J j=1 E j . ( 60
)
4 Properties of Generic Uncertainty Generic uncertainty U (p) is the length of the major semi-axis of the uncertainty ellipsoid defined by the generic covariance (C(p) in Section 2) and a probability p ∈]0, 1[. Let X 2 3 (p) be the quantile function of the X 2 distribution with 3 d.o.f. The ellipsoid is

{x ∈ R 3 , (x -p) T C -1 (p)(x -p) ≤ X 2 3 (p)}. (61)
We use notation e(p) for the smallest singular value of C -1 (p) and obtain

U (p) = X 2 3 (p) e(p) . ( 62
)
We also define reliability [START_REF] Lhuillier | Toward automatic 3d modeling of scenes using a generic camera model[END_REF]

R(p) = U (p) min i∈{1,••• ,I} ||p -o i || (63)
The topic of this Section is the study of asymptotic properties of U (p) and R(p). By "asymptotic", we mean that p should be far enough from the ray origins {o i }.

These properties rely on the main result in Section 4.1, which expresses e(p) as a function of a point in the convex hull CH of {o i }. Section 4.2 provides the link between apical angles and R, and confirms that R(p) increases linearly as the distance between p and the ray origins {o i }. Thus, U (p) increases quadratically.

Link Between e(p) and the Convex

Hull of {o i } First, we need a Lemma whose assumptions are satisfied if point p is far enough from {o i }. The proof is very technical and may be skipped by the reader.

Lemma 7 Let p ∈ R 3 such that ∃d 0 ∈ S 2 , ∀i ∈ {1..I}, d i = p -o i ||p -o i || , d T 0 d i ≥ √ 3 2 . (64) Let function G be G : d ∈ S 2 → I i=1 w i (d T i d) 2 with w i > 0. ( 65
)
There is a maximizer

d of G such that ∀i, d T i d ≥ 0.
Proof Since (a, b) is the minimal path length between a and b on S 2 , we can use properties (a, b) ∈

[0, π], π = (a, b) + (-a, b), (a, c) ≤ (a, b) + (b, c) and (a, b) = (-a, -b). Furthermore, d T 0 d i ≥ √ 3 2 and (d 0 , d i ) ≤ π 6 are equivalent. Let d ∈ S 2 such that π 3 < (d 0 , d). We have (d 0 , d) ≤ (d 0 , d i ) + (d i , d) ⇒ π 6 = π 3 - π 6 < (d 0 , d) -(d 0 , d i ) ≤ (d i , d). (66) Let d ∈ S 2 such that π 3 < (-d 0 , d). We have (-d 0 , d) ≤ (-d 0 , -d i ) + (-d i , d) ⇒ π 6 = π 3 - π 6 < (-d 0 , d) -(-d 0 , -d i ) ≤ (-d i , d) ⇒ (d i , d) = π -(-d i , d) < π - π 6 = 5π 6 . ( 67
)
Let d ∈ S 2 . Thanks to Eqs. 66 and 67, we have

π 3 < (d 0 , d) and π 3 < (-d 0 , d) ⇒ π 6 < (d i , d) < 5π 6 ⇒ |d T i d| < √ 3 2 ⇒ G(d) = I i=1 w i (d T i d) 2 < 3 2 I i=1 w i . (68) However, d T 0 d i ≥ √ 3
2 and Eq. 68 imply

G(d 0 ) = I i=1 w i (d T i d 0 ) 2 ≥ 3 2 I i=1 w i > G(d). (69)
Now, we see that all d ∈ S 2 such that π 3 < (d 0 , d) and π 3 < (-d 0 , d) are not maximizers of G. Let d be a maximizer of G (d is a singular vector of

I i=1 w i d i d T i ). We have (d 0 , d) ≤ π 3 or (-d 0 , d) ≤ π 3 . If (-d 0 , d) ≤ π 3 , -d is also a maximizer of G such that (d 0 , -d) = (-d 0 , d) ≤ π 3 . Thus, there is always a maximizer d of G such that (d 0 , d) ≤ π 3 . Last, (d, d i ) ≤ (d, d 0 ) + (d 0 , d i ) ≤ π 3 + π 6 = π 2 . ( 70
) This d is a maximizer of G such that ∀i, d T i d ≥ 0.
Here is the core of the asymptotic properties: the smallest singular value e(p) of C -1 (p) is a function of a point in the convex hull of the ray origins {o i }.

Theorem 1 Let p ∈ R 3 such that ∃d 0 ∈ S 2 , ∀i ∈ {1..I}, d i = p -o i ||p -o i || , d T 0 d i ≥ √ 3 2 . (71)
Let e(p) be the smallest singular value of C -1 (p) and

e(m, p) = 1 σ 2 α I i=1 1 ||p -o i || 2 (1 -(d T i p -m ||p -m|| ) 2 ). ( 72
(m -p) T ||m -p|| C -1 (p) m -p ||m -p|| = min m∈R 3 \{p} 1 σ 2 α I i=1 1 ||p -o i || 2 (1 -(d T i p -m ||p -m|| ) 2 ) = min m∈R 3 \{p} H(m) with H(m) = e(m, p). (74)
We also apply Lemma 7 to function

G : d ∈ S 2 → I i=1 (d T i d) 2 ||p -o i || 2 : (75) there is a maximizer d G of G such that ∀i, d T G d i ≥ 0. Since the assertions • m is a minimizer of H • p-m ||p-m|| is a maximizer of G,
are equivalent, we can choose a minimizer m of H using p-m ||p-m|| = d G . Furthermore, we have

0 ≤ d T G d i = f i (m) with f i (x) = d T i p -x ||p -x|| . (76)
There is at least one f i (m) > 0 since the maximal value of G is not 0. We define the (no null) vector

d(m) = I i=1 f i (m)d i ||p -o i || 2 . ( 77
)
The Jacobian of H is

J H (m) = 1 σ 2 α I i=1 -2f i (m)J fi (m) ||p -o i || 2 (78) with the Jacobian of f i J fi (m) = d T i (p -m) ||m -p|| 3 (p -m) T - 1 ||m -p|| d T i . ( 79
)
Since m is a minimizer of H, Eqs. 78 and 79 imply

0 = - 1 2 σ 2 α ||m -p||J T H (m) = I i=1 f i (m) ||p -o i || 2 ( d T i (p -m) ||m -p|| 2 (p -m) -d i ) = p -m ||p -m|| I i=1 f 2 i (m) ||p -o i || 2 - I i=1 f i (m)d i ||p -o i || 2 . ( 80
)
Thanks to Eqs. 80 and 77 and p

-m = ||p -m||d G , ∃λ > 0, ||p -m||d G = p -m = λd(m). (81)
Furthermore, d(m) is a positive linear combination of d i (thanks to Eqs. 77 and 76) and we obtain

∃λ i ≥ 0, d G = I i=1 λ i (p -o i ). ( 82
) Point o = p + -1 I i=1 λ i d G = 1 I i=1 λ i I i=1 λ i o i . (83
Last, Eq. 71 implies that p / ∈ CH and the proof is finished.

Asymptotic Properties

First, we show in Eqs. 88 and 90 the relation between R and apical angles (p

-o i , p -o) or (p -o 1 , p -o 2 ). According to Theorem 1, there is o ∈ CH such that e(p) = 1 σ 2 α I i=1 1 ||p -o i || 2 (1 -(d T i p -o ||p -o|| ) 2 ). ( 85
)
Let

β i = (p -o i , p -o). Since p is far from CH, sin β i ≈ β i and ||p -o i || ≈ ||p -o||. Thus, e(p) = 1 σ 2 α I i=1 sin 2 β i ||p -o i || 2 ≈ 1 σ 2 α ||p -o|| 2 I i=1 β 2 i . (86) Since U (p) = X 2 3 (p) e(p) and R(p) ≈ U (p) ||p -o|| , (87)
we see that

R(p) ≈ σ α X 2 3 (p) I i=1 (p -o i , p -o) 2 . ( 88
a 2 + b 2 σ 2 α ||p -o|| 2 ≈ 1 σ 2 α ||p -o|| 2 2(p -o 1 , p -o 2 ) 2 2 2 ≈ (p -o 1 , p -o 2 ) 2 2σ 2 α ||p -o|| 2 (89)
and we see that

R(p) ≈ σ α 2X 2 3 (p) (p -o 1 , p -o 2 ) . ( 90
)
Last, we show how U and R increase if p goes far from the ray origins o i . The double area of triangle opo i has two expressions

||(o -o i ) ∧ (p -o)|| = ||p -o i ||.||p -o|| sin β i . (91) Since sin β i ≈ β i and ||p -o i || ≈ ||p -o||, we obtain (p -o i , p -o) ≈ 1 ||p -o|| ||(o -o i ) ∧ p -o ||p -o|| ||. (92)
Eqs. 88 and 92 imply that R(p) increases linearly as the distance between p and the ray origins. Thus, U (p) increases quadratically.

Local 3D Models from Images

This Section explains how to obtain a local 3D model from a few images (I images). Local model reconstructs scene parts which are visible in a reference image.

The calibration function of the camera is known and maps image pixels to rays. It acts as a look-up table, as required in the generic context. A ray is a half line defined by its origin and direction. Thanks to the knowledge of camera pose by structure-frommotion (Section 1.4), origin o and direction d (||d|| = 1) of rays are known in the world coordinate system.

First, the reference image is segmented by a 2D mesh using gradient edges and color information (Section 5.1). Second, 3D points are reconstructed by dense stereo and ray intersection (Section 5.2). Third, geometric and reliability tests are defined from generic covariance (Sections 5.3 and 5.4). Fourth, 2D triangles are back-projected in 3D to fit the reconstructed points by tacking into account depth discontinuities, hole filling, 3D point uncertainty and reliability (Section 5.5). Last, Section 5.6 discusses the extension of these methods for non-central cameras.

2D Mesh

Two standard assumptions are used in this step. First, the scene surface should be smooth enough to be approximated by a list of triangles in 3D. Second, the occluding contours and the tangent discontinuities of surfaces are projected at gradient edges or color discontinuities in images.

The 2D mesh in the 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, uniform sampling of the field of view, and good aspect ratio for triangles. A compromise is obtained with the method described in Appendix B.

The 2D mesh has two kinds of edges: constrained edges (image contours or image borders) and unconstrained edges (Delaunay edges).

Dense Stereo and Ray Intersection

In this step, dense multi-view correspondences are calculated and reconstructed.

We apply standard pair-wise stereo method after local rectifications (Section 1.1). In our context, the local rectifications are defined by mappings from catadioptric images into faces of virtual cubes. Then the quasi-dense propagation method [START_REF] Lhuillier | Match propagation for image-based modeling and rendering[END_REF] is applied between two parallel faces of two cubes. The resulting epipolar curves are conjugate and parallel lines, except for the faces which contains the epipole: the epipolar lines intersect the epipole at the face center. Virtual cubes are preferred to virtual cylinders and catadioptric (donut) images for the reasons given in Section 1.4. In practice, cube faces are slightly extended since 3D points may be projected in two cube faces which are not parallel.

The stereo method is defined for I images as follows. We consider each catadioptric image pair (ref , sec) with ref the reference image and sec a secondary image. First, two-view dense stereo is applied for a cube of ref and a cube of sec such that faces are pair-wise parallel. Second, the stereo results are combined in the original catadioptric image ref . For each pixel of ref , the corresponding points in all secondary images sec are obtained from the matching between cube faces and the mappings between cubes and catadioptric images. The calibration function provides ray origins o i and directions d i of matched points in images i ∈ {1, • • • I}. Third, 3D points are estimated by the ray intersection method in Section 3.3. If the final value of cost function E is greater than a threshold (or if one of the I rays is not available), we can legitimately doubt the matching quality and no 3D point is retained. Small gaps (pixels of ref without 3D points) are filled with 3D points by interpolation.

Geometric Tests

Here we use the generic covariance C(p) to define several tests, which are systematically used by mesh operations for 3D modeling. Eq. 8 provides the expression of C(p) using point p, ray origins {o} and scale σ α . Section 3.3 describes a method to estimate σ α .

Let Π be the plane n ⊤ x + d = 0. The Mahalanobis point-to-point and point-to-plane [START_REF] Schindler | On robust regression in photogrammetric point clouds[END_REF] 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 . ( 93
)
The point-to-point 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). Reminder: X 2 3 (p) is defined in Eq. 61. The point-to-plane neighborhood test T (p 1 , Π) is true if d 2 (p 1 , Π) ≤ X 2
3 (p). The coplanarity 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 {p i }.

Reliability Test

A point p reconstructed from rays (o i , d i ), i ∈ {1 • • • I} may be so inaccurate that the 3D model should not contain it. Such an example occurs if camera locations are collinear: p is inaccurate and should be rejected if it is too close to the line supporting the o i .

At first glance, we can decide that p is "reliable" enough for 3D modeling if uncertainty U (p) (Eq. 62) is less than a threshold U max . However, a reliabilitybased decision is preferred: the reliability test

T (p) is true if R(p) ≤ R max with R(p) = U (p) min i ||p -o i || (94)
and R max a threshold. Now we give advantages of reliability over uncertainty for the decision.

Since the camera is central, 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. Thus, uncertainty thresholding should have a threshold U max which is proportional to the scene scale to obtain a decision which does not depend on the scale. A first reason to use Eq. 94 is its scale independence.

Furthermore, Eq. 94 allows reliable points for 3D modeling of the scene to have greater uncertainties if they are a long distance from ray origins o i and smaller uncertainties if they are close. More precisely, the permitted maximal uncertainty is proportional to the distance between point p and ray origins o i . Thus, close foreground and far background of the scene are modelized at different uncertainties. This is better than the uncertainty thresholding U (p) < U max which rejects too much background.

We note that the reliability test in Eq. 94 has two properties: (1) points in the neighborhood of the line supporting the o i (if any) are unreliable and (2) the set of reliable points is bounded. These properties result from Eqs. 88 and 92.

2.5D Mesh

Now the geometric and reliability tests (Sections 5.3 and 5.4) are used to generate a 2.5D mesh by backprojection of the 2D mesh in the reference image using the dense cloud of reconstructed points. Sections 5.1 and Section 5.2 explain how to obtain the 2D mesh and the cloud, respectively.

The principle of the method is the following. First, the 2.5D mesh is initialized as a list of fully disconnected triangles in 3D. Then, this mesh is refined by alternating operations "Triangle Connection", "Hole Filling", "Triangle Removal", "Triangle Damping" and "Mesh Refinement". Last, unreliable triangles are rejected.

At any step, the 2.5D mesh in 3D is a backprojection of the 2D mesh in the reference image: each triangle t 2d of the 2D mesh corresponds to a triangle t 3d of the 2.5D mesh with vertices v i ∈ R 3 . The t 3d vertices are parametrized 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.

In Section 5.5, the index i is used for vertex numbering, not for image numbering as in all previous Sections. Thus ∀i, o i = o since the camera is central.

Mesh Initialization A RANSAC procedure is applied to each triangle t 2d . First, all 3D points reconstructed from pixels inside t 2d are collected in a list L t 2d . Second, planes are calculated for random samples of 3 points 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, Π)} (95)
with X 2 3 (p) and d(p, Π) introduced in Eqs. 62 and 93. 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 three 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 connected 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 made 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 5.3 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 the 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 coplanarity 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.

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 a coplanarity condition and if the corresponding {t 2d a , t 2d b } in 2D have a common edge. Section 5.3 defines the coplanarity condition by T ({v i }) with {v i } the list of all triangle vertices of the k-group.

Triangle Removal A smooth surface of the scene 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, 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 foliage. 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 of ray which goes across t 3d center. In this case, the angle θ between t 3d normal n and d 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 reset by cos(θ 0 )d + sin(θ 0 ) d || d|| with d = n -(n ⊤ d)d
. "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. They 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, since the latter includes potential surface discontinuities at constrained edges too easily in the hole. Thus, 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. We would like a well defined interpolation and a low risk of depth interpolation between foreground and background. Thus we request that the hole border is coplanar using coplanarity condition T ({v i }) defined in Section 5.3. We also request enough 3D information at the hole border using 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 follows. 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. We set depth constraints 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 . Operations "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 the 3D point cloud 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 piecewise smooth surface on the 2.5D mesh.

The cost function e 3d ({z i }) is defined by

t∈T E 2 t + λ {t1,t2}∈Ed 1 2 (|t 1 | + |t 2 |)(n t1 -n t2 ) 2 (96)
with T the list of 2D mesh triangles which have triangles 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, Ed 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. 95. 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. Let z n i be the i-th depth at iteration n of the descent method. At iteration n + 1, we choose

z n+1 i ∈ {z n i -δ i (z n i ), z n i , z n i + δ i (z n i )} which minimizes the partial function z i → e 3d (z i ).
Generic uncertainty is used to scale the increment δ i by δ i (z) = ǫU (o i + zd i ) with ǫ = 0.02. 5. Remove triangles with unreliable vertex (Eq. 94).

Algorithm Summary

Step 2 connects triangles with strong conditions before step 3. Once step 3 has removed (or damped) improbable and unconnected triangles, step 4 connects triangles with weaker conditions.

From Central to 100% Generic Camera

The method in Section 5 has many limitations.

The first limitation is due to the use of a 2D mesh in a reference image (Section 5.1). The 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). Thus we should 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 should be included in the 2D mesh border. Here is a simple example of discontinuity: the line between two composite images in the generic image of a stereo-rig.

The second limitation is due to the generic covariance definition used in the geometric and reliability tests (Sections 5.3 and 5.4). The definition of C(p) requires the ray origins corresponding to p. If the camera is (approximated by) a central camera or if p is reconstructed by ray intersection, the ray origins are known. They are unknown in other (non-central) cases, unless we apply the projection functions to p (but this is not a generic method). More investigations are needed to estimate efficiently ray origins in a generic camera framework for non-central cameras.

From Local to Global 3D Models

A local 3D model (Section 5) is not adequate for a complex scene since it is view-centered. Section 6 explains how to obtain a global 3D model of the scene from a list of local models reconstructed along the sequence. Typical lists are obtained by reconstruction of one local model for each I-tuple of consecutive still images (or video key-frames) of the sequence. Global model reconstruction from local models involves view point selection, redundancy reduction, merging into topological manifold, mesh simplification, texture merging and packing. Our work only focuses on view point selection and redundancy reduction using generic covariance. Other topics are outside the paper scope.

View Point Selection

We formalize the view point selection problem as follows. Point p is reconstructed in local model l 0 and we would like to know if l 0 is one of the available local models which reconstruct p with the best accuracies. If so, p is retained in the global model. Note that p may be reconstructed in a large number of local models (in a typical list of local models) at very different accuracies, especially if the camera has a wide field of view.

Let U l (p) be the generic uncertainty (Eq. 62) of p using one ray origin o i for each image of local model l. The list of reconstructed local models is L. Local model l 0 is one of the best local models to reconstruct p if

U r l0 (p) ≤ 1 + ǫ with U r l0 (p) = U l0 (p) min l∈L U l (p) (97)
and threshold ǫ ≥ 0. In practice, ǫ > 0 is useful to reduce the lack of triangles due to matching failures (false negatives). At this point, several remarks arise. First, the view point selection defined by Eq. 97 requires the calculation of ray origins o i for p and for all local models in L. As discussed in Section 5.6, this problem is not yet solved in the non-central case. This is not a problem in our case where the camera is (approximated by) a central camera: the o i are the locations of camera center. Second, we note that Eq. 97 does not depend on probability p and scale σ α involved in the generic uncertainty definition (Eqs. 62 and 5). Indeed, changing p or σ is a multiplication of all U l (p) by the same value.

A triangle of l 0 is retained in the global model if (at least) one of its three vertices p satisfies Eq. 97. Thus, the time complexity of the view point selection method is |L| 2 |V | with |L| the number of local models and |V | the number of 3D vertices in a local model (|V | is assumed to be constant).

Note that this definition of view point selection does not depend on p visibility in views of the sequence (nor does the generic uncertainty definition). In fact, the use of U l (p) for view point selection has no sense if p can not be reconstructed by l due to lack of visibility. Thus Eq. 97 is improved by taking visibility into account as follows: we reset U l (p) = +∞ if p is not in the view field of (at least) one image of l. We may also consider the global surface to be reconstructed as a possible occluder for p in each image of l, but this problem is not integrated in this paper.

Redundancy Reduction

The method in Section 6.1 generates a global model from selected parts of local models which have the smallest uncertainties. Although the result is ready for visualization, a redundancy reduction method is useful to decrease the number of triangles.

Triangles with the largest uncertainties are progressively removed if they are overlapped (up to their uncertainty) by other triangles of the global model. A region decreasing-like method is used. At each step, we focus on the triangle t with the largest generic uncertainty which is on the border of one of the meshes currently retained in the global model. The uncertainty volume of t is defined as follows. The t vertices are v i = o i + z i d i , i ∈ {1, 2, 3} with depths z i and observations rays (o i , d i ) in the reference image of the local model l which reconstructs t. The volume is a truncated triangular cone with face f + :

v i +U l (v i )d i and face f -: v i -U l (v i )d i .
Then we test if t is overlapped by other triangles: (1) all triangles which intersect the volume are collected in a list (2) the volume is sub-sampled into segments connecting f + and f -(3) t is overlapped if all segments intersect triangle(s) in the list. If t is overlapped, it is removed from the global model. The process stops when there is no overlapped triangle in mesh borders. In practice, step (1) is accelerated using hierarchical bounding boxes and test eliminations.

Experiments

First, Sections 7.1 and 7.2 illustrate properties of generic covariance (Sections 3 and 4). Second, experiments are provided for our catadioptric cameras. The accuracy of the local 3D model reconstruction method (Section 5) is evaluated on a synthetic scene in Section 7.3. Section 7.4 illustrates the view point selection method (Section 6) to generate the global 3D model. Last, the overall system is experimented in Section 7.5 on real image sequences.

Properties of Generic Uncertainty and Reliability

Figure 1 shows generic uncertainty U and reliability R using values X 3 2 (0.9) = 6.25 and σ α = 0.001.

Black edges are the main axes of uncertainty ellipsoids centered at certain points p (their length is 2U (p)).

Every point p has a gray level depending on the interval where R(p) lies: [0, In the first case (on the left), we have two ray origins o 1 and o 2 . Both U (p) and R(p) are defined everywhere (except on the line defined by o 1 and o 2 ). Due to the symmetry of the problem, U and R are the same for any plane in 3D containing o 1 and o 2 . Furthermore, U and R increase in two cases: (1) if p goes toward a line defined by o 1 and o 2 or (2) if p goes far from {o 1 , o 2 }. We also see at the bottom that our generic reliability is similar to the angle-based reliability used in [START_REF] Doubek | Reliable 3d reconstruction from a few catadioptric images[END_REF][START_REF] Torii | Measuring camera translation by the dominant apical angle[END_REF]: curves implicitly defined by constant R(p) are very similar to circles defined by constant apical angle (o 1p, o 2p). In the second case (in the middle), we add a ray origin o 3 in the middle of o 1 and o 2 . The result is unexpected: there is no improvement (i.e. U or R decrease) by adding o 3 . In fact, the results are nearly the same.

In the third case (on the right), we slightly move o 3 toward the bottom. As expected, the improvement is noticeable in the neighborhood of the line defined by o 1 and o 2 . In these two last cases, our R definition is naturally derived from U for any numbers of views. This is not the case for angle-based reliability [START_REF] Doubek | Reliable 3d reconstruction from a few catadioptric images[END_REF], which is only defined for two views.

Last, a numerical test is provided for the asymptotic relation (Eq. 90) between reliability R and apical angle (o 1p, o 2p) in two cases: two ray origins o 1 and o 2 , and three ray origins with o 3 = 1 2 (o 1 +o 2 ). Function

a(p) = (o 1 -p, o 2 -p) R(p) σ α 2X 3 2 (p) (98)
should have a small standard deviation σ a and a mean ā close to 1 if p is far enough from ray origins.

Values σ a and ā are estimated from samples p in a ring centered on o 3 with large radius r l = 10||o 1 -o 2 || and small radius r s = 1 2 ||o 1o 2 || (the ring is in a plane which contains the ray origins). We obtain (ā, σ a ) = (1.06, 0.15) in the two ray origin case and (ā, σ a ) = (1.05, 0.11) in the three ray origin case.

Larger r s and r l improve the results, as expected.

Errors in Unit Sphere, Image and 3D

We have seen that the generic covariance definition is based on S 2 error (Section 2). Furthermore, the projection function p of the camera propagates this error to image error such that the uncertainty ellipses of image error are distortion ellipses of p (Section 3.1). Figure 2 shows these ellipses for our equiangular catadioptric camera (on the right) and two other cases: a perspective projection into a cube face (on the left) and an equirectangular projection (in the middle). The equirectangular projection maps 3D point to its spherical coordinates (ϕ, θ) and is often used for panoramic imaging (Figure 2 only shows a half part of the view field). We see that the propagated image error depends on p and is not "standard" (isotropic and uniform in the whole image), especially for cameras with wide view fields. We deduce from Figure 2 that a local rectification from the catadioptric image (on the right) into cube face (on the left) provides a more standard image error. A yet more standard image error is obtained if we replace the perspective projection (left of Figure 2) by the equirectangular projection (right of Figure 3).

Figure 3 illustrates the virtual cubes and the epipolar geometry involved in the dense stereo step of our reconstruction method (Section 5.2). The equirectangular projection is used for the local rectifications into cube faces. According to Section 3.2, this choice makes the generic covariance similar to the virtual covariance.

Accuracy of Local 3D Models for a Synthetic Scene

Now, quantitative results are given for a local 3D model reconstructed from catadioptric images: scene accuracy (discrepancy between scene reconstruction and its ground truth) and calibration accuracy (discrepancy between the calibration and its ground truth). Summary First, synthetic images are generated using ground truth (non-central) calibration. Second, structure-from-motion [START_REF] Lhuillier | Automatic scene structure and camera motion using a catadioptric system[END_REF] (SfM) and the 3D modeling method described in Section 5 are applied using many central calibrations defined by radial distortion functions. Third, the scene reconstructions are registered with ground truth and quantitative evaluations are made. Details on non-central calibration and registration are shown later in Appendix A. The scene reconstructions are degraded for many reasons: image noise, approximate calibration, small baseline (due to distant or collinear points with the camera motion) and low textured areas.

Experiment Choice The catadioptric camera moves inside a textured cube to experiment the reconstruction in the whole field of view. A representative range of baselines is obtained with the following ground truth: the [0, 5] 3 cube and camera locations defined by

o i = 1 1 + i 5 1 ⊤ , i ∈ {0, 1, 2}
(numbers in meters). Camera orientations (rotations) are slightly perturbed around I 3 to simulate hand-held camera. The catadioptric image is a ring with large radius of 1128 pixels and is projected into cube faces such that each face has about 440 × 440 pixels (as in real experiments).

Scene Accuracy A standard definition of error err(p) for a reconstructed point p is the signed distance of p to the ground truth surface. In the context of a camera rotated around a small object, it is a natural choice to tolerate a uniform error for all parts of the object [START_REF] Seitz | A comparison and evaluation of multi-view stereo reconstruction algorithms[END_REF]. In our context of a camera moving in a scene including close foreground and far background, it is more adequate to define an error tolerance which increases with point depth. Thus, our accuracy measure a 0.9 is defined by the 90% fractile of |err(p)| ||p-o1|| for all triangle vertices p of the 3D model. In other words, |err(p)| ≤ a 0.9 ||po 1 || is true for 90% of vertices.

Experiments Figure 4 shows reconstruction results of the synthetic cube using three central calibrations "Best", "Shift" and "Sfm". These calibrations are defined by radial distortion functions β(r), which Figure 4: Reconstruction results of a synthetic cube using three calibrations: "Best" (row 1), "Shift" (row 2), and "Sfm" (row 3). Calibrations are defined by radial distortion functions which map radius r of image point to angle β between ray and camera axis. From left to right: mapping r → β(r) superimposed with its ground truth, error graph, triangle vertices and faces of ground truth cube projected onto 3 planes (more details in the text). All numbers are in centimeters, except in the first column. Row 4 shows many texturemapped views of the local 3D model reconstructed with "Sfm" and one image of the sequence. Each local 3D model is constructed from 3 camera locations, which are shown in the middle column. map normalized radius r of image point to angle β between ray and camera axis (small and large circles are r = 0 and r = 1, respectively).

The first column shows many β(r) graphs. Camera poses are estimated using Structure-frommotion [START_REF] Lhuillier | Automatic scene structure and camera motion using a catadioptric system[END_REF], and local 3D models are obtained using method of Section 5. "Best" is a very accurate approximation of radial distortion defined by the true (non-central) calibration. "Shift" is "Best" shifted by 2 degrees. "Sfm" is estimated by structurefrom-motion [START_REF] Lhuillier | Automatic scene structure and camera motion using a catadioptric system[END_REF]. The standard deviation of β (with ground truth) of calibrations "Best", "Shift" and "Sfm" are 0.077, 1.98 and 0.65 degrees, respectively.

The second column of Figure 4 shows the corresponding error graphs defined by joint distribution of depth x = ||po 1 || and error y = err(p) for all triangle vertices p of the 3D models. 90% of dots (x, y) are between the two lines y = ±a 0.9 x with a best 0.9 = 0.0085, a shif t 0.9 = 0.033 and a sf m 0.9 = 0.015. In examples "Shift" and "Sfm", β(r) inaccuracy produces a majority of vertices inside the ground truth cube where error err(p) is negative. The β(r) inaccuracy also increases a 0.9 . The cloud dispersion of "Sfm" is slightly larger than that of "Best".

The other columns of Figure 4 project the local 3D models into planes (xy), (xz) and (zy), which are parallels to faces of the ground truth cube. In the "Shift" case, calibration inaccuracy produces camera pose distortion, cube distortions and a majority of vertices inside the ground truth cube. These problems are less visible for "Sfm".

The row on the bottom shows texture-mapped views and triangle orientations of the reconstructed model "Sfm". A large neighborhood of a cube corner (on the left) and a small part of the red ground (on the right) are not reconstructed since they are not in view field.

Local 3D Model Selection for Catadioptric Camera

The global 3D model of a scene is obtained by view point selection of local 3D models (Section 6.1): a point p reconstructed in local model l 0 is retained in the global model if the relative uncertainty U r l0 (p) is less than a threshold (Eq. 97). Section 7.4 illustrates how this method works for a typical list of local 3D models (one local model is reconstructed for each triplet of consecutive images of the sequence).

Figure 5 shows values of U r l0 for a catadioptric camera pointing toward the sky and moving on the horizontal ground with almost collinear camera locations. In the first case (columns 3 and 4 on the right), the local model is not at the end of the whole sequence. We note that a kind of planar slice of the 3D space contains small values of U r l0 (p), with U r l0 (p) = 1 at the central component. The planar slice goes across the middle camera, its thickness increases with the distance to the middle camera, and it is connected to both ends of the whole camera sequence. It should be remembered that visibility is used in the view point selection as follows: we reset U l (p) = +∞ if p is not in the view field (of at least one image) of local model l. In the catadioptric case, we have U l (p) = +∞ if p is in a blind cone of l shown on the top of Figure 5. This has several consequences: (1) U r l0 is not defined outside the view field of l 0 , (2) U r l0 is not continuous at boundaries of view fields of l with l = l 0 , and (3) the reconstruction of ground surface is allowed. Figure 5 shows these 3 consequences in column 4:

(1) white-gray checker-board, (2) gray level discontinuities and (3) low values of U r l0 in the immediate neighborhood of the camera trajectory on the ground (if the shifted plane of the figure is the ground surface). In the second case (columns 1 and 2 on the left), the local model is at the end of the whole sequence. The planar slice is replaced by a large section of a half space, and the three notes on visibility are still correct.

Figure 6 shows results of view point selection (Eq. 97 with ǫ = 0.1) combined with unreliable point rejection (Eq. 94 with R max = 0.1, σ α = 0.001 and X 2 3 (0.9) = 6.25). On the left, gray levels encode the minimal uncertainty available to reconstruct scene point p with all local 3D models of the typical list L of a 11-view sequence. On the right, gray levels encode the number of local models available to reconstruct scene point p thanks to view point selection. This number is the modeling redundancy at p. We note that redundancy increases if p leaves the camera trajectory, and redundancy is limited thanks to unreliable point rejection.

Real Examples

Now, the overall reconstruction system is experimented on still images sequences taken by equiangular catadioptric cameras, which are hand-held and mounted on a monopod (such a camera choice is discussed in Section 1.3). We obtain 3264 × 2448 JPEG images with the following setup: adequate mirrors [1] mounted with the Nikon Coolpix 8700 using adapter rings (Figure 7). Both image dimensions are divided by 2 to accelerate structure-from-motion (SfM), mesh estimations, and to facilitate texture storage of the VRML models. Only dense stereo benefits from original image sizes using the local rectifications into cube faces (Section 5.2). Church Sequence We take 208 images during a complete walk around a church using the 360 One VR (model 3) mirror, using a supplementary adapter ring which adds 2 cm between mirror and perspective camera. The trajectory length is about (25 ± 5cm) × 208 = 52 ± 10m (the exact step lengths between consecutive images are unknown). The radii of large and small circles of the catadioptric images are 1128 and 232 pixels.

SfM is the first step of the method (Section 1.4). Figure 8 shows top views of the SfM results obtained with slight modifications of the method [START_REF] Lhuillier | Automatic scene structure and camera motion using a catadioptric system[END_REF]: the number of inlier points is bounded by 500 in each image to accelerate hierarchical bundle adjustment and (optional) loop closure, then a last global bundle adjustment is applied without inlier bounds. The unclosed reconstruction and its drift are shown in the top left corner; its closed version is shown in the top-right corner. The final result is shown at the bottom with some images. It has 76033 3D points and 477744 points in images satisfying the geometry. The final RMS error is 0.74 pixels. A 2D point is considered as an outlier if its reprojection error is greater than 2 pixels. With our setup, the estimated view field angles are α up = 41.5, α down = 141.7 degrees (the angles given by the mirror manufacturer are α up = 40, α down = 140 degrees).

The second step is the calculation of a list of local 3D models (Section 5). Here we use typical list: one local model for each triple of consecutive images of the sequence. Figure 9 shows the reference image, the depth map and the 2D mesh of a local model. The reference image is between the two others. Once the catadioptric images are projected into cube faces, quasi-dense propagation [START_REF] Lhuillier | Match propagation for image-based modeling and rendering[END_REF] is used and benefits by the large number of seed matches provided by SfM inliers. The depth map is given in the original reference image. The 2D mesh has 22894 triangles and the 2.5D mesh is generated using X 2 3 (0.9) = 6.25 and σ α = 0.00082 radians. This value of σ α is estimated from the dense cloud reconstruction (using Eq. 60) of all local models of the list. Figure 9 also shows views of the local model without and with rejection of unreliable triangles (using Eq. 94 and R max = 0.05). The local model with rejection has 10510 triangles. A small part of the ground 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.

Once the 208 local 3D models of the closed sequence are reconstructed (with the same value of σ α ), the global 3D model is obtained by view point selection and redundancy reduction (Section 6). First, triangles of local models are selected using the reliability test: all vertices p of a triangle should satisfy Eq. 94 with R max = 0.04. We obtain 2249675 triangles. Second, view point selection is applied: any triangle of any local model is retained in the global model if it has at least one vertex p such that Eq. 97 is true with ǫ = 0.1. Only 31% of triangles are retained (681561 triangles). The (final) global model has 368814 triangles after redundancy reduction. Figure 10 shows many views of this global model of the church and its neighborhood. The reader can match the top view of the model (top of Figure 10) with the top view of the SfM result in Figure 8. Now, the difficulties of this scene and the consequences on the global 3D model are shown in detail. There are four trees in the immediate neighborhood of the church and many others which are more distant. Since trees are important scene components, we choose to not remove triangles which are unconnected with others to modelize the foliage as clouds of textured triangles. Thus, the "Triangle Damping" mesh operation is preferred and used instead of the "Triangle Removal" operation (Section 5.5). The drawback is the presence of certain isolated and false triangles in other scene components. Furthermore, the scene has many low textured areas: streets, cars, and sky. The consequences are the presence of a number of holes and inaccurate reconstructions in the streets or distant cars. Reconstruction is almost impossible for the blue car which is very close to the camera trajectory (see picture in the top-left corner of Figure 8). The sky is globally unmatched thanks to the confidence measure used by quasi-dense propagation [START_REF] Lhuillier | Match propagation for image-based modeling and rendering[END_REF], but minor parts of the sky occur in reconstructed foliage and at the top of some buildings. Last, the images are taken in the early morning. Two consequences are: (1) the shadow of the author occurs in several images and (2) the first and last images of the closed sequence have different textures. The former Old Town Sequence This sequence has 354 images and is acquired using the 0-360 mirror. The radii of large and small circles are 1140 and 204 pixels. The trajectory length is about (35±5cm)×353 = 122 ± 17m.

SfM estimates 149792 3D points with 819087 2D inliers and RMS error equal to 0.75 pixels. These inlier numbers are larger than those in [START_REF] Lhuillier | Toward flexible 3d modeling using a catadioptric camera[END_REF] thanks to a slight modification of the SfM method. Figure 11 shows a top view of the SfM result and some images of the sequence. The estimated view field angles are α up = 34.5, α down = 152.9 degrees (the angles given by the mirror manufacturer are α up = 37.5, α down = 152.5). Then, the 352 local models of the typical list are estimated using X 2 3 (0.9) = 6.25 and σ α = 0.00089 radians. The numbers of triangles after unreliable triangle rejection (R max = 0.03), view point selection (ǫ = 0.1) and redundancy reduction are 4991183, 1786919 and 1184620, respectively. Figure 12 shows many views of the final global model. The joint video is a scene walkthrough in the scene.

We note that the reconstruction of building tops is inaccurate since it involves several difficulties. First, building tops are not visible by the closest local model but by others which are more distant (small baseline). Second, there is a lack of calibration accuracy (angle α up ) in the neighborhood of the large circle where building tops are projected . Third, false matches are not rejected by epipolar constraint since the camera motion is parallel to the building borders. In this context, unreliable triangle rejection is useful to avoid these parts in the global model. Other difficulties are low textured areas (buildings, ground, cars), specular reflections on windows, and a major illumination change during image acquisition.

We also note that a few parts of the scene are not retained in the global model by view point selection, although they are not in a blind cone of the local model which provides the smallest uncertainty. These parts sometimes include walls with plane normal parallel to the camera trajectory at the closest camera location. This is due to our current and limited use of visibility in the view point selection (end of Section 6.1).

Conclusion

This paper is a new step towards the automatic 3D modeling of scenes using a generic method, which can be applied for any kind of camera. First, we introduce a generic covariance for central cameras.

Proof is given for its definition and several properties: asymptotic behaviors, links with the apical angles, links between the generic covariance and the covariance which results from the ray intersection problem defined by the sum of squared reprojection errors in pixels. Second, we present a 3D modeling method from images which systematically uses the generic covariance. More precisely, generic covariance is used to weight minimized scores involved in triangle estimation, to decide connections between triangles, to fill holes, to reject unreliable triangles, to define view point selection and to reduce model redundancy.

However, these contributions do not go far enough for non-central camera. The calculation of generic covariance at point p requires the ray origins corresponding to p. These ray origins are known if the camera is (approximated by) a central camera, but a generic method is still needed to obtain ray origins for a non-central camera. Other limitations are due to the use of 2D meshes in images: the camera should not be too exotic to back-project connected image points (2D triangles) to connected points of the scene surface. Furthermore, the current implementation of 2D mesh initialization is only done for our cameras. All these limitations are topics for future research.

Experiments are another contribution of the paper. Our context is useful for applications and has rarely been addressed until now: the automatic 3D modeling of scenes using catadioptric cameras. Once structure-from-motion has estimated camera calibration and poses from image sequence, the 3D modeling method based on generic covariance is applied to generate a global 3D model of the scene. Experiments illustrate generic covariance properties. They include accuracy estimation of our catadioptric setup on syn- Many improvements are possible: a better use of visibility in view point selection, applying merging methods like [START_REF] Curless | A volumetric method for building complex models from range images[END_REF][START_REF] Soucy | A general surface approach to the integration of a set of range views[END_REF] after view point selection. It is important to remember that view point selection is a key issue for 3D modeling, but it can not replace (and it is not a competing method of) such merging methods. Lists of local models with different baselines would also be useful to increase accuracy. Last, several steps of the reconstruction method may be improved (image matching, texture merging) and other steps may be added (mesh simplification, image gain corrections).

Appendix A: Data for Synthetic Experiments

This appendix explains how to obtain a realistic noncentral calibration of a catadioptric camera. This is useful for the synthetic experiments in Section 7.3.

Here we use the 0-360 mirror [1] mounted with the Nikon Coolpix 8700. Non-central calibration is required to project the synthetic scene in images by taking into account ray reflexion on the mirror. It is defined by the mirror profile and the matrix of the perspective camera in the coordinate system of the mirror [START_REF] Lhuillier | Automatic scene structure and camera motion using a catadioptric system[END_REF].

First, the mirror profile has to be estimated since it is not provided by the mirror manufacturer. Profile z(r) is used by the mirror cylindrical parametrization f (r, θ) = r cos(θ) r sin(θ) z(r) ⊤ (99) and is obtained as follows: (1) select a room with lighting reduced to a single bulb on the ceiling (2) put a Cartesian graph paper on the ground such that the light rays are orthogonal to the paper (3) put the mirror on the paper such that the mirror symmetry axis is horizontal and parallel to one of the main directions of the Cartesian graph paper (4) mark several points on the paper at the border of the mirror shadow and (5) fit polynomial z(r) from these points. We obtain z(r) = 0.0186 + 0.06188r + 0.10812r 2 + 0.0312154r 3 with 0 ≤ r ≤ 3.9 cm. (100) Second, the matrix of the perspective camera is estimated. A value of the focal length f p is available from the EXIF data in the JPEG image returned by the perspective camera: f p = 7094 pixels. Furthermore, a typical radius of the large circle in catadioptric images is r up = 1128 pixels. Following Section 4.5 of [START_REF] Lhuillier | Automatic scene structure and camera motion using a catadioptric system[END_REF], we obtain the z-coordinate of the perspective camera center on the mirror symmetry axis using the Thales relation 3.9 z(3.9)+zp = rup fp . Now, we obtain the perspective camera matrix K p R ⊤ p I 3 -t p in the mirror coordinate system: intrinsic parameters K p = diag(f p , f p , 1), rotation R p = I 3×3 and perspective center t p = 0 0 -z p ⊤ with z p = 20.8 cm.

Third, we check the view field defined by these values. We obtain α up = 37.4 using ray reflection at the large circle. In practice, the ratio between radii of small and large circles in the catadioptric image is equal to 0.18. This ratio and ray reflection at the small circle imply α down = 153.0. These angle values are similar to those of the mirror manufacturer (α up = 37.5 and α down = 152.5 degrees).

We also define the 3D registration which maps the estimated scene reconstruction to scene ground truth. This is useful since the definition of reconstruction accuracy in Section 7.3 depends on it. This registration is a similarity transformation defined as follows. Its Euclidean part maps the local basis of the 2 nd pose of the reconstruction (central camera model) to the local basis of the mirror at the 2 nd pose of the ground truth (non-central camera model). Its scale part is the ratio of trajectory lengths between ground truth and reconstruction bases. 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 e 2d ({p v }) is defined by

p∈t∈T ||c p - p ′ ∈t c p ′ |t| || 2 + λ v∈V || v ′ ∈Nv p v -p v ′ || 2 (101)
and 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 vertex locations {p v } are the parametrization of e 2d , and e 2d is minimized using a descent method. 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 priority to detected gradient edges over the minimization of color variance, which may sometimes be contradictory.

)

  There is o(p) in the convex hull CH of {o i } such that e(p) = min m∈CH e(m, p) = e(o(p), p). (73) Proof Thanks to the definition of e(p) and the expression of C(p) in Eq. 8, e(p) = min m∈R 3 \{p}

  ) is both in the convex hull CH of {o i } and in the line defined by direction d G and point p. Point o is a minimizer of H, as all other points in this line (except p). Thanks to this definition of o and Eq. 74, we obtain e(p) = min m∈CH\{p} e(m, p) = e(o(p), p).

  ) Similarly, if there are two ray origins, Theorem 1 implies e(p) = min o∈[o1,o2] e(o, p) ≈ min |a|+|b|=(p-o1,p-o2)

7 20 π) 4 .

 74 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 = 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 two ray origins (left), three collinear ray origins (middle) and three non collinear ray origins (right). Ray origins 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. On the left, black curves are circles defined by constant apical angles.

Figure 2 :Figure 3 :

 23 Figure 2: Distortion ellipses for the perspective projection into a face cube (left), the equirectangular projection (middle) and our equiangular catadioptric camera (right). The circular cone aperture 2ǫ of all distortion ellipses is π 25 .

Figure 5 :

 5 Figure 5: Top: three locations (black points) of a catadioptric camera moving in horizontal plane (view fields are bounded by cones). Middle: mapping p → U r l0 (p) on this plane and a shifted plane with l 0 two 3-view local 3D models of a 15-view sequence. Bottom: local zooms of figures in the middle. The camera locations of local model l 0 are black points; others locations are white. Values [1, 1.5] of U r l0 are mapped to colors white-gray; larger values are dark gray; undefined U r l0 (p) area is white-gray checkerboard.

Figure 6 :

 6 Figure 6: Left: gray levels encode p → min l∈L U l (p) = U l(p) (p) in the plane where the camera moves. Areas with the same best local model l(p) are bordered by black lines. Right: gray levels encode the number of local models accepted by view point selection (darkest gray: 1 local model, white: 8 local models). Black pixels on the image borders are points rejected by reliability condition.

Figure 7 :

 7 Figure 7: 360 One VR (left) and 0-360 (right) mirrors mounted with the Nikon Coolpix 8700.

Figure 8 :

 8 Figure 8: Top: unclosed and closed structure-from-motion results with reduced number of 3D points. Bottom: closed structure-from-motion with all 3D points and some images of the Church.

Figure 9 :

 9 Figure 9: Top: reference image and depth map of a 3-view local model. Middle left: local view of the 2D mesh (constrained and unconstrained Delaunay edges are red and blue, respectively). Middle right: local model without reliability test. Bottom: local model with reliability test (triangle orientations on the right).

Figure 10 :

 10 Figure 10: From top to bottom: top view and height map of church global model, local view (texture and orientation of triangles).

Figure 11 :

 11 Figure 11: Top view of the structure-from-motion result with some images of the old town.

Figure 12 :

 12 Figure 12: Left: top view and height map of the old town global model. Right: local views (texture and orientation of triangles).

Appendix B: 2D Mesh

This appendix is a summary of the 2D mesh method used in Section 5.1. There are three steps.

Mesh Initialization First, a Delaunay triangulation is initialized in the reference image such that the solid angles of all triangles are roughly the same. This step is defined as follows for our catadioptric cameras (the general case has been left out for future work).

We generate a checkerboard such that (1) each cell is quadrilateral and has two triangles (2) the mesh resolution is defined by a mean length of cell edges equal to 8 pixels (3) catadioptric image borders enforce constrained edges on the Delaunay and enforce the global shape of the checkerboard. At this point, cell rows of the checkerboard are concentric rings and cell columns are radial sections. A modification of the current mesh is useful to increase the uniformity of solid angles of triangles: we remove one vertex out of two at image border neighborhoods to increase the solid angles of triangles in these areas.

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 gradient edges are integrated in a best first order. A contour is a list of connected pixels which have maximum local image gradient. The contour 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 these 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" operations (move vertices to minimize a global cost combining color variance in triangles and mesh smoothness) and "discrete" operations (flip edges and merge vertices to improve aspect ratio of triangles).

The continuous operation is useful for many reasons. First, a few 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