Maxime Lhuillier

Toward Flexible 3D Modeling using a Catadioptric Camera

Fully automatic 3D modeling from a catadioptric image sequence 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, and the majority of them require calibration of the camera. This paper presents a method which deals with hundreds of images, and does not require precise calibration knowledge. In this context, the same 3D point of the scene may be visible and reconstructed in a large number of images at very different accuracies. So the main part of this paper concerns the selection of reconstructed points, a problem largely ignored in previous works. Summaries of the structure from motion and dense stereo steps are also given. Experiments include the 3D model reconstruction of indoor and outdoor scenes, and a walkthrough in a city.

Introduction

Producing photo-realistic 3D models for walkthroughs in a complex scene given an image sequence 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 horizontal plane, when the viewer moves along the ground. This suggests a wide field of view for the given images, for which many kinds of cameras are possible [START_REF] Daniilidis | The page of omnidirectional vision[END_REF]: catadioptric cameras, fish-eyes, or systems of multi-cameras pointing in many directions. Since we would like to capture any scene (indoor and outdoor) where a pedestrian can go, the hardware involved should be hand-held/head-held and not cumbersome. A catadioptric camera is a good candidate for all these constraints, and it has been adopted in this work.

The main drawback of this choice is the low resolution compared with a standard (perspective) camera for a given field of view. We would like to compensate for this problem using still image sequence taken with an equiangular catadioptric camera. Still images are preferred to video images thanks to 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 in the whole image. These two choices mainly have two consequences. First, a still image sequence requires some effort and patience: the user should alternate a step forward and a (non blurred) shot. Second, an equiangular catadioptric camera is not a central camera. A non-central model complicates all 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 a such setup, and that a central approximation for the camera is sufficient in many cases.

The presented approach is fully automatic given an image sequence acquired by the camera moving in the scene:

(1) estimate the full geometry using a Structure from Motion (SfM) method, and (2) build a 3D model using multiview dense stereo and stereo fusion. By SfM, we mean the automatic, robust and optimal estimation of the camera motion (all extrinsic and some intrinsic parameters) and the scene structure (a sparse set of points in 3D) from images. The 3D model is a list of textured triangles in 3D which approximate the visible part of the scene, and is visualized using standard tools. Although this approach is now "standard" for a monocular, hand-held and perspective camera (e.g. [START_REF] Nister | Automatic Dense Reconstruction from Uncalibrated Video Sequence[END_REF][START_REF] Pollefeys | Visual modeling with a handheld camera[END_REF][START_REF] Lhuillier | A quasi-dense approach to surface reconstruction from uncalibrated images[END_REF]), no such 3D reconstruction results have been obtained before with a catadioptric camera.

Previous Works

During the past decade many researchers have worked on the 3D model reconstruction given an image sequence. A lot of work has been done in different contexts involving active or passive sensors, aerial or pedestrian view points, general or specific camera motions (e.g. turntable), general or specific objects to modelize (e.g. faces), manual or automatic methods. A complete survey is outside the scope of this paper, and we only focus on the most related work.

Previous authors [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] reproject their original images onto virtual cylinder (panoramic) images to apply dense stereo with the epipolar constraint. The first attempt [START_REF] Kang | 3-d scene data recovery using omnidirectional multibaseline stereo[END_REF] merges images taken by a rotated perspective camera, and the others [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] reproject calibrated (central) catadioptric images. Point tracking and the 8-point algorithm [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] are used by [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] to estimate the geometry of pair of successive cylinder images. Then, the relative scales between geometry pairs are recovered by measuring the distance between camera poses or odometry [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], and the global geometry of the sequence is obtained. Finally, these two approaches use multi-baseline stereo methods inspired by [START_REF] Okutomi | A multiple-baseline stereo[END_REF]. The third work [START_REF] Fleck | Omnidirectional 3d modeling on a mobile robot using graph cuts[END_REF] uses active methods to recover the sequence geometry, and a graph-cut method [START_REF] Kolmogorov | Computing visual correspondances with occlusions using graph-cuts[END_REF] followed by a lot of postprocessing for dense stereo. A simple manual method is also proposed by [START_REF] Doubek | Reliable 3d reconstruction from a few catadioptric images[END_REF], but the important point here is the introduction of an image pair selection scheme to reconstruct a point with the best "reliability". This is a key issue (neglected by other authors) for 3D modeling from catadioptric images, since the same 3D point may be visible and reconstructed in a great number of images at very different resolutions and baselines.

As mentioned before, SfM [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] is the first problem to solve. Although the principles are well known, optimal (including global bundle adjustment) and robust SfM systems are not so common if the only given data are a long image sequence acquired by a general catadioptric camera and some knowledge about the calibration. This contrasts with the situation regarding perspective cameras. The most advanced research into this subject have been [START_REF] Strelow | Extending shape-from-motion estimation to noncentral omnidirectional camera[END_REF] and [START_REF] Micusik | Structure from motion with wide circular field of view cameras[END_REF]. A bundle adjustment (BA) is applied once at the very end of the geometry estimation process given an accurate camera calibration [START_REF] Strelow | Extending shape-from-motion estimation to noncentral omnidirectional camera[END_REF], although it is recognized that a hierarchical scheme for BA [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF] greatly improves the accuracy and robustness simultaneously. The uncalibrated approach [START_REF] Micusik | Structure from motion with wide circular field of view cameras[END_REF] is a two step camera modeling method: (1) recover the sequence geometry using an approximate central model for the camera and (2) upgrade the sequence geometry using a non-central model enforcing the mirror knowledge.

Contributions

Section 2 presents a central camera model and summarizes the SfM method. The differences with [START_REF] Micusik | Structure from motion with wide circular field of view cameras[END_REF] are: long sequences thanks to a hierarchical BA, general radial function with some rough parameter knowledge (instead of the knowledge of model type with unknown intrinsic parameters). Section 3 presents a method to obtain a "local" 3D model of the scene given three images of the sequence, using dense stereo. We don't use cylinder images like [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] in the dense stereo step to facilitate the 3D reconstruction of the scene ground. Section 4 describes a generalization of the image pair selection scheme to reconstruct a point: we replace the heuristic two-view criterion [START_REF] Doubek | Reliable 3d reconstruction from a few catadioptric images[END_REF] (using angle between rays) by a general multi-view criterion (using uncertainty ellipsoids). The most accurate parts of all local models are selected by this criterion, and they are merged into the global 3D model of the scene. This subject was neglected by other authors, although the same 3D point may be estimated with very different accuracies (in our catadioptric context) depending on the images selected for reconstruction. Note that selection is a preprocessing step for 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], it is not a concurrent method. Experiments in Section 5 include both indoor and 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] Doubek | Reliable 3d reconstruction from a few catadioptric images[END_REF][START_REF] Fleck | Omnidirectional 3d modeling on a mobile robot using graph cuts[END_REF]).

Geometry Estimation (Summary)

In this Section, the automatic SfM approach is briefly summarized (more technical details in [START_REF] Lhuillier | Automatic structure and motion using a catadioptric camera[END_REF]). The following assumptions are required: (1) a surface-of-revolution mirror whose lower and upper circular cross sections are visible (2) a perspective camera with zero skew and aspect ratio set to 1 (3) mirror and perspective camera with the same symmetry axis (4) constant calibration and (5) approximate knowledge of the two field of view angles.

Central Camera Model

The camera model is defined by its orientation R (a rotation), the center t ∈ 3 ( is the real numbers), both expressed in the world coordinate system, and a central projection p : 3 \ {0} → 2 . Using the same notation X for a finite 3D point and its world coordinates, the direction of the ray from t to X in the camera coordinate system is given by d = R (Xt). The image projection of X is p(d). The model has a symmetry around the z-axis of the camera coordinate system: the omnidirectional image is between two concentric circles of radii r up and r down , and there is a positive and decreasing function r such that p(x, y, z) = r(α(x, y, z))

x √ x 2 +y 2 y √ x 2 +y 2 with α(x, y, z) the angle between the z-axis and the ray direction d = x y z . Let α up and α down be the two angles which define the field of view. We have α up ≤ α ≤ α down and r down ≤ r(α) ≤ r up . An exact equiangular camera is obtained if r is a linear function.

Calibration Initialization

Thanks to the assumptions, the projections of the mirror circular cross sections are concentric circles. First these circles are detected and estimated in each catadioptric image using RANSAC and Levenberg-Marquardt methods applied on regularly polygonized contours (in this paper, the detection stability is improved assuming that the circles are fixed in the whole sequence). Then, the initial r(α) is defined by the linear function such that r(α up ) = r up and r(α down ) = r down . The image circle radii r down , r up are obtained in the first step, and the field of view angles α down , α up are given by the mirror manufacturer. These angles are not exactly known since they depend on the relative position between the mirror and the pinhole camera.

Geometry Initialization

Harris points are detected and matched for each pair of consecutive images in the sequence using ZNCC correlation, without any epipolar constraint. The corresponding ray directions are also obtained from the calibration initialization. Then, the essential matrices for these pairs are estimated by RANSAC (using the 7-point algorithm [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]) and refined by Levenberg-Marquardt. 3D points are also reconstructed for each pair. Many of these points are tracked in three images, and are used to initialize the relative 3D scales between two consecutive image pairs. Finally, the full sequence geometry is obtained by many bundle adjustments (BA) applied in a hierarchical framework to merge all partial geometries [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. These BA simultaneously refine the parameters of 3D points and camera poses by minimizing the sum of squared reprojection errors.

Geometry Refinement

Once the full geometry is obtained for the approximate function r(α) defined above, r(α) is redefined as a cubic polynomial whose the 4 coefficients should be estimated. An additional BA is applied to estimate the 4 + 6c + 3p parameters of the sequence (c is the number of cameras, p is the number of 3D points) and increase the numbers of 3D and 2D inliers.

Local 3D Models

Assume that the geometry of the catadioptric image sequence is estimated with the method described in Section 2. Since a central camera model is used, the epipolar constraint and reprojections to virtual surfaces (like cylinder or cube) are easy. This is not the case with a non-central model.

A local model is reconstructed from three images: one reference image ref , one secondary image sec 1 before and one secondary image sec 2 after the reference image in the sequence. The three images do not have to be consecutive to deal with many baselines and reconstruct parts of the scene at various depths. Obviously, large baselines increase the accuracy for distant parts, but small baselines are also necessary for the closest parts (e.g. ground) to increase the common field of view and simplify the matching problem.

Previous authors [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] reproject their images onto virtual cylinder images to apply dense stereo with the epipolar constraint. The resulting epipolar curves are sinusoids [START_REF] Mcmillan | Plenoptic modelling: an imagebased rendering system[END_REF]. The main disadvantage of cylindrical images is the high image distortion for top and bottom parts, corresponding to the sky and the ground. This increases the difficulty of the matching problem for these scene components. For this reason, we prefer to reproject a catadioptric image onto the 6 faces of a virtual cube and apply a classical two-view dense stereo (for convenience [START_REF] Lhuillier | Match propagation for imagebased modeling and rendering[END_REF]) to two parallel faces of two cubes. The resulting epipolar curves become parallel lines, except for the faces which contain the epipoles: the epipolar lines intersect the epipole at the face center.

Thus, the method for a catadioptric image triple

(ref ,

Global 3D Model

Let L be the list of all local models that we have computed with the methods in Sections 2 and 3. Some parts of these models should be selected and merged into a global 3D model. First, Section 4.1 defines the virtual uncertainty U l (P ) as a function of any 3D points P and any local models l of L. The virtual uncertainty is the usual uncertainty if P is reconstructed by l, and the virtual uncertainty extends it for other local models. Second, the local 3D model selection for reconstruction is presented in Section 4.2. Given a point P reconstructed by a local model l 0 , we wish to know if l 0 is one of the local models in L with the smallest virtual uncertainties for P . If this is not the case, P should not be retained in the global model since a better model is available for P reconstruction. Even if l 0 minimizes the virtual uncertainty for P , the quality of P may be too bad for a global model (e.g. if P and all camera centers of l 0 are collinear points). For this reason, some "reliability" conditions for P are described in Section 4.3. Last, Section 4.4 describes how to use the selection criterion efficiently and obtain a global model from all local models of L.

Virtual Uncertainty Definition

The virtual uncertainty U l (P ) for a 3D point P and a local model l is defined as follows. Let J l be the Jacobian of the function 3 → 2k projecting P on all images i 1 • • • i k from which l is reconstructed. We note that C l (P ) = σ 2 (J l (P )J l (P )) -1 is the covariance matrix associated to the reconstruction problem of P from im-ages i 1 • • • i k , assuming independent and identical Gaussian noise of σ pixels for reprojection errors and zero uncertainties for cameras. An estimate of σ is provided by the reprojection errors of ray intersections (see Section 3). We define the virtual uncertainty U l (P ) as the length of the major semi-axis of the uncertainty ellipsoid defined by covariance C l (P ) and a probability p.

This definition of U l (P ) assumes that P is "visible" in all images of l. If this is not the case, we consider that l can not reconstruct P accurately and choose U l (P ) = +∞. The visibility definition depends on the knowledge we have of the global surface to be reconstructed. Obviously, a visible point P should be in the view fields of all images of l. Furthermore, the surface normal at P defines a half-space which should contain each camera location of l. We may also consider the global surface to be reconstructed itself as a possible occluder for P in each image of l, but this problem is not integrated in this paper.

The virtual uncertainty extends the "reliability" proposed by [START_REF] Doubek | Reliable 3d reconstruction from a few catadioptric images[END_REF]: only the case k = 2 was considered and the reliability was defined by π 2 -arccos(|d 1 .d 2 |), with d 1 and d 2 the directions of the rays which go across P and camera centers i 1 and i 2 respectively. Although intuitive (best/worst reliabilities are for two perpendicular/parallel rays), the reliability does not depend on the distance between camera centers. Uncertainty does.

Local 3D Model Selection for Reconstruction

The selection criterion is defined using the virtual uncertainty. Assume that a point P is reconstructed by a local model l 0 . We would like to know if l 0 is one of the best local models in L to reconstruct P .

One can estimate U l (P ), ∀l ∈ L and sort them in increasing order. If U l0 (P ) is ranked at top n, l 0 is one of the best local models. A threshold n > 1 is useful to increase the density of points retained in the global model, and also to tolerate certain matching failures (false negatives). However, the virtual uncertainties at top n may have different magnitude orders. We avoid this drawback by estimating the relative uncertainty U r l0 (P ) = U l 0 (P ) min l∈L U l (P ) and retain l 0 as one of the best models for P if U r l0 (P ) ≤ 1 + . As has been previously stated, a threshold > 0 is useful to increase the point density in the global model.

We note that these criterion selections are independent of the choice of probability p and noise σ used to define the virtual uncertainties: changing p or σ is just a multiplication of all U l (P ) by a same scale factor.

Reliability Conditions

Assume that a point P is reconstructed by a local model l 0 . Point P may be considered as unreliable and rejected if it is too far away from the camera centers of l 0 , or if these centers and P are collinear points. These cases are possible even if l 0 is known to be one of the best local models. So, a reliability condition is useful to decide whether P should be included in the global model or not.

At first glance, the reliability condition might be "P is reliable if U l0 (P ) < U 0 " with U 0 a threshold. However, the rendering of the resulting global model will suffer from two drawbacks in our visualization purpose (walkthrough in the scene). First, points at the foreground will have greater reprojection errors than points in the background. Second, too many points in the background will be discarded since point uncertainty increases with depth. So, a second definition for the reliability condition might be "P is reliable if U l0 (P ) < U 0 d γ (P, l 0 )" with γ > 0 and d the mean distance between P and the camera centers of l 0 . Now, the two drawbacks are reduced, but this definition is heuristic and depends on thresholds p (a probability defined in Section 4.1), U 0 and γ.

A third definition is given by [START_REF] Doubek | Reliable 3d reconstruction from a few catadioptric images[END_REF]: P is reliable if there is a camera pair (i, j) such that π 2 -arccos(|d i .d j |) < π 2θ 0 with θ 0 a threshold and d i the direction of ray which goes across P and the i-th camera center of l 0 . This definition requires only one threshold θ 0 (a lower bound for angle between two camera rays). It deals with both cases "P is too far away" and "P and camera centers of l 0 are collinear points".

From Local to Global Model

A straightforward use of the selection criterion for the global model generation has a time complexity proportional to (#L) 2 (#P ) with #L being the number of local models in L and #P the (average) number of 3D points in a local model. This complexity may be high since #L is at least proportional to the image sequence length (sometimes many hundreds) and #P has the same magnitude order as the number of pixels in a catadioptric image (hundred of thousands).

The time calculation is reduced by subdividing each local model into small patches, and by applying the selection test (and the reliability test) only one time for each patch given a representative 3D point. Once a global set of patches is selected from all local models, usual merging/fusion methods [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] may be used to reduce the redundancy of overlapped patches in space. In this paper, each patch is a square (assembled in rings of catadioptric images) and its representative point has the median depth of the patch. Although the reconstruction system does not yet include a merging method, the recovered 3D models are convincing as shown in the next Section. Two textured triangles are used for each patch.

Experiments

After the presentation of the experimental context in Section 5.1, synthetic experiments on local 3D model selection are presented in Section 5.2. Sections 5.3 and 5.4 give and discuss results obtained with the 3D scene reconstruction system for real sequences.

Context

The user moves along a trajectory on the ground with the catadioptric camera mounted on a monopod, alternating a step forward and a shot. A roughly equiangular catadioptric camera is used (the "0-360" mirror with the Nikon Coolpix 8700, shown in Figure 1). This is a non-central camera, such that the symmetry axes of camera and mirror are assumed to be the same. The mirror manufacturer gives the view field angles α up = 37.5 and α down = 152.5 degrees. These angles are not exactly known since they depend on the relative position between the mirror and the pinhole camera. Image dimensions are 1632 × 1224 pixels.

Local 3D Model Selection for Reconstruction

This part presents the experiment of the local 3D model selection (Section 4) with synthetic camera motions. First, the uncertainty map U l (P ) is shown on the left of Figure 2 for two local models. As expected, U l (P ) increases if point P goes near the line where all camera centers of l lies (if any) or if P goes away from camera centers of l. The first case (local model with 3 collinear view points) often occurs if the given camera trajectory is a smooth curve on the ground. Since this local model reconstructs many points close the given trajectory with a great deal of uncertainty, these points will appear very noisy in a rendered image even in the favorable context of a rendering view point in the neighborhood of the given camera trajectory. Selecting points in the global model with low uncertainty is a poor solution since distant points will also be rejected. High uncertainty is more acceptable for distant points than close points if the rendering view points is in the neighborhood of the given trajectory. The second local model (with 3 non-collinear view points) is a more favorable case since P is never aligned with all camera centers of l.

The relative uncertainty U r l (P ) is shown on the right of Figure 2 for 3-view local models in a 7-view sequence with nearly collinear camera centers. As described in Section 4, a point P reconstructed by l should have a small U r l (P ) to be accepted in the global model: U r l (P ) should be in [1, 1 + ]. In the first case, 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 l (P ), with U r l (P ) = 1 at the central component. The planar slice goes across the middle camera of the local model, its thickness increases with the distance to the middle camera, and it is connected to both ends of the whole camera sequence. In the second case, the local model is at the end of the whole sequence and the planar slice is replaced by a large section of a half space. In both cases, we see that the local 3D model selection (i.e. thresholding U r l (P )) is not sufficient to decide whether the point of a local model should be accepted in the global model.

Figure 3 gives results obtained by combining the local 3D model selection and the angle-based reliability condition defined by θ 0 (see Section 4.3). In all cases, black pixels outside the gray-white region are points which do not satisfy the reliability condition with θ 0 = 10 o . On the left, two circles are the set of points P such that the angle between the two camera centers are equal to θ 0 . In the middle and on the right, we consider all 3-view (consecutive) local models of an 11-view sequence. The space partition defined by the best local model is given in the middle. On the right we see the number of local models satisfying the local model selection = 0.1 with θ 0 = 10 o for a given P (i.e. the number of possible reconstructions for P ). We note that this number increases if P goes out of the camera trajectory. The second step is the calculation of all local 3D models (see Section 3). A local model is build for each triple of consecutive images of the sequence, and is defined by a list of about 10000 patches partitioning the reference image. Figure 4 patches (819767 patches) are retained in the global model.

Old Town Sequence

Figure 6 shows three views of this model. It is not difficult for the reader to match the oblique view of the full model (top of Figure 6) with the top view of the SfM result in Figure 5. Furthermore, the joint video is composed of two parts. The first part shows all original images. Panoramic images are also given to help understand the scene, but they are not used by the method. A major lighting change occurs at frame number 295. The second part is a walkthrough in the whole scene in the neighborhood of the original view points. We consider that the results are acceptable in textured areas in spite of the low resolution (about 210-260 pixels for an angle in space of π 4 ). Currently, the main problem is due to many parts of walls and streets which are not textured enough to be matched and reconstructed. A second problem is due to the patches which have regular size (about 8 × 8 pixels) and locations in images. The resulting approximation of occluding contours is sometimes very crude.

Quantitative results for the 90% point uncertainties are given in Figure 7 for the global model. In 99, 5 percent of cases, the uncertainty increases from 0.6 cm to 40 cm while point depth grows from 1 m to 9 m.

Small Room Sequence

The Small Room sequence is composed of 18 indoor images, such that the step length between two consecutive view points is about 15 cm (exact values are unknown). SfM results are: 2974 reconstructed points, 15096 image points satisfying the geometry with RMS=0.73 pixels, We note that a single local model at the middle of the sequence is not sufficient to reconstruct all objects, and more local models are welcome to reduce the uncertainties.

Conclusion

A fully automatic 3D modeling method from a catadioptric image sequence is proposed. First, still images are acquired by a roughly equiangular camera. Second, the sequence geometry (including some intrinsic parameters) is successfully estimated with a central camera model. Third, many local 3D models along the sequence are reconstructed with 3 views. The image reprojection on a virtual cube (instead of a cylinder) allows the use of classical dense stereo methods and facilitates the reconstruction of the ground. Last, the global 3D model is obtained by applying a local model selection: each local model is partitioned in patches, and a patch is rejected if an other local model is available to reconstruct the patch with less uncertainty. Such 3D models have never been obtained before with catadioptric cameras, because matching is always difficult in practice and the last step has been largely ignored in previous studies.

Many improvements are possible and include: better use of visibility in the local model selection, a better choice of patches near occluding contours, improving matching in low textured areas, and more investigations on the choice of local models to calculate (decrease their number, integrate different baselines for a same reference image).

Figure 1 .

 1 Figure 1. The "0-360" mirror with the Nikon Coolpix 8700.

Figure 2 .

 2 Figure 2. Uncertainty maps P → U l (P ) or U r l (P ) in the horizontal plane (smaller values are white). The camera centers of local model l are black points in this plane. From left to right: U l (P ) for 3 aligned cameras, U l (P ) for 3 non-aligned cameras, U r l (P ) for a 3-view local model in the middle of a 7-view sequence, U r l (P ) for a 3-view local model at the beginning of a 7-view sequence. The neighborhoods of camera centers are zoomed in the second row using histogram normalization. Right: the black color is used if 1.5 < U r l (P ), and camera centers are white if not in l.

Figure 3 .

 3 Figure 3. Left: the uncertainty map P → U l (P ) of a 2-view local model l. Middle and right: all local models of 3 consecutive views are considered in a 11-view sequence. Middle: the uncertainty map P → U l(P ) (P ) with l(P ) the best local model (the local model minimizing l → U l (P )). The regions of points P which have the same l(P ) are bordered by black lines. Right: the number of local models l accepted by the local model selection ( = 0.1) is encoded by a gray level (darkest gray: 1 local model, white: 9 local models). In all cases, pixels are black if the corresponding P does not satisfy the reliability condition with θ0 = 10 o .

Figure 4 .

 4 Figure 4. Reference image 211 and its depth map (white pixels have no depth).

  It requires about 52 minutes to take the 354 images of the Old Town sequence. Some images are shown in Figure 5. The trajectory length is about (35 ± 5cm) × 353 = 122 ± 17m (the exact step lengths between consecutive images are unknown). The radii of large and small circles of the catadioptric images are 570 and 102 pixels. A rectangular field of view of π 4 × π 4 has approximately 260 × 210 pixels in the horizontal and vertical directions in theses conditions. SfM is the first step of the method as described in Section 2. 59859 3D points are automatically reconstructed with 380947 points in images satisfying the geometry. A top view of the result is proposed in Figure 5. The final RMS error is 0.74 pixels, and a 2D point is considered as an outlier if the reprojection error is greater than 2 pixels. All calculations are started with an inaccurate calibration (α up = 40, α down = 140 degrees) to experiment the robustness of the SfM and the behavior of the calibration refinement. The recovered radial distortion r(α) is very close to a linear function (as expected) and the view field angles are improved: α up = 35.5, α down = 152.2 degrees.

Figure 5 .

 5 Figure 5. Some images of the Old Town sequence and a top view of the reconstruction using SfM, including 354 cameras (black squares) and 59859 points (black points).

  shows the depth map and the reference image of a local 3D model. Last, a global 3D model is obtained by selecting patches of all local models using the local model selection scheme (Section 4.2) and a reliability condition (Section 4.3). A patch of a local model is accepted if its representative 3D point (introduced in Section 4.4) satisfies (1) the angle-based reliability condition with θ 0 = 5 o and (2) the local model selection with = 0.1. Only 24% of

Figure 6 .

 6 Figure 6. Three textured-mapped views of the Old Town global model (including cameras), obtained with local 3D model selection and reliability conditions. Patch orientations (or depth map) are also drawn using gray levels.

Figure 7 .

 7 Figure 7. Quantitative results for uncertainty of global models (Old Town on the left, and Small Room on the right). Each patch has coordinates (d l (P ), U l (P )) with P the representative point of the patch, d l (P ) and U l (P ) the mean depth and the uncertainty (p = 90%) of P by the local model l which reconstructs P . The ranges (cm) are [103, 890] × [0.63, 39.9] (left) and [55.6, 212] × [0.72, 7.56] (right), and exclude 0.5% of patches.

α

  up = 37.8 and α down = 152.7 degrees (with α up = 40 and α down = 140 as initial values). Only 26% of patches (35709 patches) are retained in the global model after local model selection ( = 0.1) and reliability constraint (θ 0 = 5 o ).

  Figure 7 gives quantitative uncertainties for the global model and Figure8shows many views of the global model. The main objects are easy to recognize. As mentioned in the Old Town example, many parts of the room are not reconstructed due to the lack of texture, and the patches ignore the occluding contours.

Figure 8 .

 8 Figure 8. Four textured-mapped views of the Small Room global model, obtained with local 3D model selection and reliability conditions. Depth map and patch orientations are drawn on the top and in the middle. One image of the sequence is also shown.

  sec 1 , sec 2 ) is as follows. First, two-view dense stereo is applied for a cube of ref and a cube of sec 1 , and for a cube of ref and a cube of sec 2 . Second, the stereo results are combined in the original catadioptric image ref .For each pixel of ref , the corresponding points in sec 1 and sec 2 are obtained from the mappings between cubes and the mappings between cubes and catadioptric images. A 3D point is obtained for the current pixel of ref by intersection of three rays using Levenberg-Marquardt. If the reprojection error is greater than a threshold (or if one of the three rays is not available), we can legitimately doubt the matching quality and no 3D point is retained. Last, many gaps (pixels of ref without 3D points) are filled with 3D points by interpolation.


Bernard De Gouvello

Pierre-Antoine Versini

Julien Paupardin

Manuel Pruvost- Bouvattier

MULTI-SCALE APPROACH OF THE GREEN ROOFING POTENTIAL: APPLICATION TO THE PARIS AGGLOMERATION

Green roofs are experiencing a significant development in France since the early 2000s. In 15 years, the market has increased from 100,000 m² implemented / year to over 1.3 million. If the motivations for the implementation of green roofs are often multiple, their participation in urban storm water management issues and in urban heat islands reduction are particularly often cited (Li and Babcock, 2014).

. It is based on the combination of systematic land cover data for the three considered scales (provided by the IAU-IF: Institut d'Aménagement et d'Urbanisme d'Ile-de-France) crossed with urban development and history of the regional architecture and the use of a statistical approach based on a sample of buildings. The calculation and representation of the GRP is adapted on the studied territorial level. The results shows that the GRP differs significantly from a county to another and also within a same county.

INTRODUCTION

The development of green roofs (GR) is increasing in France for ten years. Around 100,000 m² were installed per year in the early 2000s, and now the market is around 1.3 million m² (data Adivet, 2013). GR have often been presented as offering benefits for the city in terms of aesthetic contribution, thermal regulation, biodiversity support, air quality improvement and storm water management. A significant number of studies were focused on this last topic over the past decade. They aimed to characterize and/or model the hydrological behavior of green roofs (eg: Villareal and 1-4 December 2015 2 Bengtsson, 2005; [START_REF] Berndtsson | The influence of extensive vegetated roofs on runoff water quality[END_REF][START_REF] Baraglioli | Studies on the quantity impacts of green roofs[END_REF][START_REF] Palla | Unsaturated 2D modelling of subsurface water flow in the coarse-grained porous matrix of a green roof[END_REF][START_REF] Berthier | Simulation of green roof hydrological behavior with a reservoir model[END_REF][START_REF] Gregoire | Effect of a modular extensive green roof on stormwater runoff and water quality[END_REF][START_REF] Stovin | The hydrological performance of a green roof test bed under UK climatic conditions[END_REF][START_REF] Gromaire | Incidence of extensive green roof structures on the quantity and the quality of runoff watersfirst results from an experimental test bench in Paris area[END_REF].

In a dense urban area, the proportion of roofs can represent a significant fraction of the surface. In such environments, characterized by low land availability which limits the implementation of "traditionnal" BMPs, GR may play a role, since they do not need additional space. Few works were focused on this issue [START_REF] Carter | Vegetated roofs for stormwater management at multiple spatial scales[END_REF][START_REF] Palla | Modelling storm water control operated green roofs at the urban catchment scale[END_REF]Versini et al, 2015b and2015A.). Based on of real urban study basins located in the Hauts de Seine, Versini et al. (2015b) showed that GR could have a significant effect to reduce the network overflow frequency or runoff discharge from sewer systems when GR is significantly implemented. The ability of a given territory to host a large area of GR becomes an issue that should be investigated further. Indeed, this capacity is highly variable from one territory to another because of the diversity of urban characteristics and particularly its "built heritage". This communication is intended as an original contribution to the analysis of this issue. Its purpose is to offer a simple definition of the potential for revegetation of an urban space and to present a method to estimate it. This method is defined and applied in the territory of Ile-de-France to three successive levels: municipality (city of Créteil), county (Seine-Saint-Denis) and the Parisian agglomeration (with more 400 municipalities of the Paris region). At each considered scale, some specific modes of representation and analysis are conducted.

The paper is organized in two parts. The first one presents the concepts and the adopted methodological approach. In the second one, the results obtained for each considered scale are presented through tables and maps, and discussed.

CONCEPTS AND METHODOLOGY

The potential for green roofing (GRP) of a given urban area is the cumulative area of the existing roofs, which are already vegetalized or technically able to be vegetalized in the future. The criterion adopted for considering a roof technically able to be vegetalized is its flatness. Indeed, it is assumed that: (i) current technologies allow to provide GR solution for any type of flat roof -hence with a seal -and (ii) the greening of inclined roofs is so marginal at urban scale that it can be neglected.

The systematic identification of flat roofs for an entire city is a time consuming and expensive task. Although it was conducted by the city of Paris in order to define a GR policy [START_REF] Apur | Etude sur le potentiel de végétalisation des toitures terrasses à Paris[END_REF], it is not possible to lead such a study on larger territories. Several methods of automatic recognition of roofs' flatness have also been developed, but they are not yet fully operational and remain costly.

General Principle

In this work, the methodology uses a different logic. It is inspired by the approach proposed by Emna Sellami-Kaaniche to evaluate pollutants flows from roofs on the urban scale and applied in the city of Créteil [START_REF] Sellami-Kaaniche | A new method for modelling roofing materials emissions on the city scale: Application for zinc in the City of Créteil (France)[END_REF]. The urban space is divided into 4 periods corresponding to the historical urban development of the Paris region: before 1900, 1900-1960, 1960-1982 and 1982-2012. The land use patterns corresponding to built-up areas are grouped into 13 distinct classes: individual homes, collective housing, other types of housing, secondary activity, offices, indoor sports, teaching equipment, health equipment, cultural facilities, administrative equipment, other equipment, commercial areas, and transports. This leads to a partition of up to 52 entities. For a given portion of territory, the value of an entity is provided by the total of ground surfaces of the buildings belonging to the same class and located on the area corresponding to the same historical period of urban development.

Each entity i has a unique flat roof ratio (TRi, %), that represents the proportion of terraces (flat roofs) among the roofs included in this entity. These ratios have been determined on the basis of samplings selected for each entity on the city of Créteil (considered representative of the urbanization of the Paris area). They have been completed for entities insufficiently represented in this city by some additional samplings in Paris and several municipalities from the Seine-Saint-Denis and Seine-et-Marne counties. Table 1 shows all the 52 applicable TRi ratios.

Calculation and mapping of GRP for the three scales

The method allows to calculate GRP for each of the entities located in any fraction of a given territory. Depending on the scale, it is necessary to conduct an adequate grouping of entities. The adopted methodology for the calculation and the mapping of the GRP depending on the studied scales is specified below.

At the municipal level (city of Créteil), the GRP calculation provides only a single value. Therefore, it is interesting to calculate the value of GRP for several sub-areas of the city. This can be done by dividing the communal space into sectors corresponding to an administrative reality. Creteil is divided into 5 sectors grouping 20 neighborhood councils' jurisdictions. It has been decided to represent on the same map the entities grouped in intervals of TR values and the GRP values of the five sectors. Creteil has also served as the basis for the validation of the method, since automated identification of flat roofs was conducted by the IAU-IF for this city.

When the scale is reduced, a more radical aggregation is necessary to achieve interpretable maps. For both county (40 municipalities) and agglomeration (431 municipalities) scales, the choice was to aggregate data inside each municipality.

At the County scale, two calculations have been conducted. First, the current GRP has been computed on the basis of buildings existing in 2012 (considered as the current situation). Second, a prospective calculation has been made. It is based on some already known future developments that are likely to change significantly the result. This is the case of the "Grand Paris" which will conduct to create railway and subway stations for new or existing lines in the surrounding counties of Paris as for the Seine-Saint-Denis. It has been estimated that such a new station may impact the urban development within 800 m. The contribution to GRP induced by the implementation of these stations has been re-calculated for each municipality. This calculation was made by considering the following assumptions: (i) the urban densification in new stations' areas will increase the buildings' footprint ratio and the corresponding new buildings will all have flat roofs; (ii) the footprint ratio of an affected area is likely to reach 0.30, without exceeding the double of the current situation ratio. The 0.30 value was established by analyzing the corresponding values for the 30 existing train stations (RER) of the County. It corresponds to the average value of the 6 highest values.

These two assumptions lead to the following calculation for each municipality:

𝐺𝑅𝑃 𝑃 = 𝐺𝑅𝑃 𝐶 + 𝑀𝑖𝑛(𝐶𝑆𝐴 * 𝐹𝑅𝐵𝑆𝐴, 𝑀𝑎𝑥[𝐶𝑆𝐴 * (0.30 -𝐹𝑅𝐵𝑆𝐴),0]) where: GRP P (ha) is the prospective GRP of the municipality; GRP C (ha) is the current GRP of the municipality; CSA (ha) is the cumulated surface of the Stations' Areas of the municipality; FRBSA is the Footprint Ratio of the Buildings in the Stations' Areas of the municipality.

Finally, at the agglomeration scale, GRP is recalculated for each municipality and is mapped. Table 2 summarizes the approach to compute and represent GRP for every level.

Table 2: GRP calculation and representation at every level

Scale

RESULTS AND DISCUSSION

Municipal scale: the city of Créteil

The GRP computed for the city of Créteil is about 100 ha. This represents approximately 8% of the city area (11.4 km²) and more than 53% of the cumulative buildings footprints (178 ha). Regarding these figures, green roofing can affect the appearance of the city. However, this potential is characterized by a high spatial variability. Some differences appear significant as for the center and the eastern sector, where a quite low potential is noticed. Conversely both western and southern sectors are characterized by a GRP more than three times higher. The northern sector represents an intermediate situation with a GRP of 16.6 ha. The difference is partly explained by the number of apartment buildings in both western and southern areas contrasting with a majority of individual housing in central and eastern sectors. In the southern sector, it is worth noting the high contribution to the potential provided by the Créteil Soleil shopping center located east of the lake. This terrace is currently used as a car park: the actualization of this potential therefore depends on the building of a superstructure.

Map 1: RT dissemination and GRP by sector for the city of Créteil

Créteil allows to assess the accuracy of the method in terms of spatial resolution. Table 3 provides comparative GRP calculation by using: (i) the presented method and (ii) a numerical method implemented by the IAU-IF and combining satellite and geographic data. This table shows a quite good matching between both methods at the municipality scale. It seems logical since the TR ratios were calculated on the whole city. At the sectors' level, computed GRP remain convincing with some similar values obtained for the three main areas (center and east, north, west and south).

Scale of County: the Seine-Saint-Denis

The current GRP of the Seine-Saint-Denis County is about 13 km². It represents more than 5.5% of County total area (236 km²). According to the prospective scenario described in the previous section, this potential is likely to increase, and might reach 21.5 km², i.e. almost 10% of the territory. This result shows that the "Grand Paris" project, which will lead to a significant densification of the built environment in this County, is a tremendous opportunity to noticeably increase the presence of vegetation in the city by the implementation of green roofs. The current situation reflects a significant GRP variation from one municipality to another. There is no simple explanatory factor (as a "center-periphery" effect for instance). The relative importance of GRP in the north-western part of the County is due to the concentration of both collective housing and industrial buildings characterized by some quite large (roof) areas.

The prospective map shows that the future railway stations are likely to strengthen significantly the GRP of many municipalities from north-west towards respectively north-east and south. Some municipalities located in the holes of the future rail transport network and on the limit of the county keep their situation unchanged, leading to a widening gap. It should however be noticed that these municipalities have mostly large green spaces regarding their total surface.

Paris agglomeration

Table 4 summarizes the results of GRP for all the Paris agglomeration and at each County level. 1-4 December 2015

It is noteworthy that the GRP represents 3.33% of the total surface of the agglomeration, but this proportion differs from a County to another. Three Counties have a greater proportion than 5%: Paris town (75), Hauts-de-Seine (92) and Seine-Saint-Denis (93). It is consequently important to follow a deepening analysis for these territories.

At a more detailed level, it appears that the situation varies noticeably from one municipality to another (see map 3).

Map 3: GRP of the Municipalities of Paris agglomeration

This map shows a very uneven distribution of values. Three-quarters of the municipalities of the metropolitan area (320) have a GRP less than 30 ha when only 7 of them exceed 100 ha. Moving away from the capital, one can notice that the GRP generally tends to decline. This can be explained by the gradual "de-densification" of the urban spaces. However, there are several notable exceptions to this trend. On the one hand, GRP is important near the international airports (Roissy and Orly). On the other hand, one can note of the GRP increase within the 5 "new towns" created in the late 1960s' and now included in the scope of the Paris agglomeration (Marne-la-Vallée, Cergy-Pontoise, Saint Quentin-en-Yvelines, Evry and Senart). These areas are characterized both by a local re-densification and a greater proportion buildings with flat roofs, corresponding to their construction era. indirect levers-players [START_REF] Belmeziti | Impact Potentiel de l'Utilisation de l'Eau de Pluie dans le Batiment sur les Consommations d'Eau Potable à l'Echelle Urbaine[END_REF]. These aspects have to be analyzed in order to develop public policies to move from assessment of GRP to RGRP (Reachable Green Roof Potential). The identification of territories with a high GRP is nevertheless a key-step in this process. Depending on the local operational issues of stormwater management (network overflow, or spill into the environment, that are likely to intensify with climate change), a deep analysis of RGRP in these territories has to be led. As such, the GRP calculation approach is a tool to support future developments of management at the source of rainwater.

Figure 1 :

 1 Figure 1: General principle of the method

  ha) calculated gathering GRP values of the 39 existing entities -Mapping of TR by gathering entities, and GRP by sector (5 sectors for Créteil) County (Seine-Saint-Denis)

  Maps 2 (a) and (b) represent the current GRP and the prospective GRP for all the municipalities of the Seine-Saint-Denis County. Maps 2: GRP of the municipalities of Seine-Saint-Denis County: (a) Current GRP; (b) Prospective GRP

Table 1 : Terraces' ratios by entity (TR i , %)

 1

	Classe de bâtiments	Période n°1	Période n°2	Période n°3	Période n°4
		(Avant 1900)	(1900-1960)	(1960-1982)	(1982-2012)

Table 3 : Comparative calculation of GRP (our method vs IAU-IF method) 1. Our method 2. IAU-IF method Difference GRP (ha) GRP (ha) (GRP1-GRP2)/GRP2 (%)

 3

	Centre	8.65	8.724	0.82
	Est	8.37	9.38	-10.79
	Nord	15.92	17.34	-8.16
	Ouest	29.46	35.49	-16.98
	Sud	31.11	25.38	22.57
	Total	84.87	96.31	-11.88

Table 4 : GRP of the Counties of Paris agglomeration

 4

			Total Area (km²)	Roofs area (ha)	GRP (ha)	GRP/Total
									Area (%)
	County Number	of	Total	Moyenne/	Total	Moyenne/	Total	Moyen/
	municip.			commune		commune		commune
	75	20	105.42	5.27	1073	53.66	649	32.47	6.16
	92	36	175.57	4.88	1347	37.40	938	26.04	5.34
	93	40	236.86	5.92	2025	50.63	1320	33.01	5.57
	94	47	244.83	5.21	1464	31.15	1153	24.54	4.71
	F-77*	52	420.55	8.09	1623	31.21	1155	22.21	2.75
	F-78*	89	678.59	7.62	2089	23.47	1470	16.52	2.17
	F-91*	82	562.84	6.86	2001	24.41	1446	17.63	2.57
	F-95*	65	434.81	6.69	1631	25.09	1470	22.61	3.38
	Total	431 2859.47	6.63	13253	30.75	9522	22.09	3.33
	Agglo.
		(

*: These Counties are partially included in the Paris Agglomeration)

ACKNOWLEDGMENTS

This work has been supported by the French ADEME framework program through the TERRACES project and by the County of Seine-Saint-Denis.

CONCLUSION

The storm water management within cities in a climate change context creates pressure on the availability of space to implement mitigation techniques solutions, such as alternative techniques. Because they do not occupy additional space, green roofs are a particularly useful tool in this context. However, this work has shown that the physical characteristics of urban space influence sensitively the possibility of installing green roofs.

The ability to revegetate the roofs is not limited solely to physical characteristics of urban space. It is also based on the consideration of legal, economic and the identification of direct and