Maxime Lhuillier

Surface reconstruction from a sparse point cloud by enforcing visibility consistency and topology constraints

Keywords: Surface Reconstruction, Manifold, Genus, Visibility, Sparse Features, Environment modeling

Surface reconstruction under

Introduction

The automatic surface reconstruction of an environment from an image sequence is still an active research topic. Here we present a method that generates a surface using topology constraints and given a sparse cloud of 3D points reconstructed from the images. First Sec. 1.1 describes these constraints and their importance, since they are under-explored in Computer Vision. Then Secs. 1.2 and 1.3 summarize previous works that reconstruct a surface from such a cloud or topology constraints, respectively. Last Sec. 1.4 presents our contributions.

Why enforcing topology constraints ?

There are two reasons to enforce constraints on the computed surface. First this regularizes the problem and helps to deal with noisy and lacking input points. Second topological constraints are needed in downstream processing and applications most of the time. There are several topology constraints. sphere (with dotted equator) meets g = 0. Middle: simple torus meets g = 1. Right: double torus meets g = 2.

Manifoldness

A surface is manifold if every of its point has a small neighborhood in the surface that is mapped to R 2 by a homeomorphism (a bijective and continuous function whose inverse is continuous). If a surface point does not meet this condition, the point is singular. Now we briefly remind arguments of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF]. Method implicitly assumes that surface is manifold if it estimates normal or curvature. This holds not only for surface reconstruction methods that regularize using smoothing ( [START_REF] Hernandez Esteban | Silhouette and stereo fusion of 3d object modeling[END_REF][START_REF] Vu | High accuracy and visibility-consistent dense multiview stereo[END_REF]), but also for applications like rendering ( [START_REF] Botsch | Polygon mesh processing[END_REF]). The manifold constraint can be seen as the weakest smoothing that can be applied to a surface, even when the density of input points is too low to estimate reliable normal or curvature. According to experiments of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF], the manifold constraint is necessary in our context: 25% of the vertices at the surface are singular if we do not enforce it.

Low genus

Assume that we have a manifold surface in R 3 . Its genus g is the number of surface handles, e.g. g = 0 for a sphere and g = 1 for a torus as shown in Fig. 1.

We say that there is topological noise ( [START_REF] Wood | Removing excess topology from isosurfaces[END_REF]) if the estimated surface has a g excess compared to the true surface, i.e. if the estimate has spurious handle(s). Both lack of points and point noise in the input are reasons of topological noise. Enforcing a low genus during surface computation can provide a strong regularization ( [START_REF] Zou | Topology-constrained surface reconstruction from cross-sections[END_REF]), but most methods are unable to do this.

The topological noise is usually removed by a postprocessing. This is useful for downstream processing such as smoothing, refinement and texture mapping ( [START_REF] Wood | Removing excess topology from isosurfaces[END_REF][START_REF] Botsch | Polygon mesh processing[END_REF]), which require low topological noise. An example is the dense stereo refinement: the surface evolves in 3D such that it minimizes a photo-consistency function without changing the genus ( [START_REF] Hernandez Esteban | Silhouette and stereo fusion of 3d object modeling[END_REF][START_REF] Vu | High accuracy and visibility-consistent dense multiview stereo[END_REF][START_REF] Romanoni | Automatic 3d reconstruction of manifold meshes via delaunay triangulation and mesh sweeping[END_REF]). The lower the topological noise, the better the expected convergence. (Think about a torus that tries to approximate a sphere by minimizing a photo-consistency score.)

Last we provide theoretical remarks, which may motivate to consider genus in Computer Vision methods that reconstruct surface. First g and Gaussian curvature K meet S K = 4π(cg) if S is a manifold surface with c connected components ( [START_REF] Dyn | Optimizing 3d triangulations using discrete curvature analysis[END_REF]). This is consistent with intuition: we can reduce the (negative Gaussian) curvature of the surface by removing handles if c g, and this looks like curvature-based regularization. Second there are assumptions ( [START_REF] Amenta | Surface reconstruction by voronoi filtering[END_REF]) on the input point cloud and the true scene surface such that a 3D Delaunay triangulation based surface reconstruction method can provide the good g.

Connectedness

Similarly as in the genus case, both noise and lack of points can generate spurious connected components of the surface. This is also topological noise. Thus connectedness can be used as constraint ([30]).

Previous sparse methods

The sparse methods estimate a surface from a 3D point cloud, which is sparse since it is computed from interest points and/or contour points having uneven and low distributions in the images. Nevertheless the points have decent confidence thanks to standard pipeline including point selections, robust and optimal reconstruction by using RANSAC and non-linear least squares.

We shortly remind advantages of the sparse methods. First they are interesting for both time and space complexities, especially for obtaining compact models of large scale environments. This is useful for computations with limited hardware resource, or for applications that need high scalability and do not need the level of details provided by dense stereo. Second they can initialize dense stereo methods to improve accuracy and obtain more detailed reconstructions, if the experimental conditions are favorable to get enough texture in the images.

Most sparse methods use a 3D Delaunay triangulation of the input points. This discretizes the space by a set T of tetrahedra. Furthermore, every point has visibility information: a set of line segments called rays linking the point and every camera location that contributed to reconstruct the point. A bipartition of T is obtained: the free-space is the set of the tetrahedra intersected by at least one ray, the free-space complement is the matter. The set of triangles separating free-space and matter is a candidate surface ( [START_REF] Lovi | Incremental free-space carving for real-time 3d reconstruction[END_REF][START_REF] Ling | Building maps for autonomous navigation using sparse visual slam features[END_REF]), but it is neither manifold nor robust to bad points. This surface can be converted into a manifold surface by vertex duplications ( [START_REF] Gueziec | Cutting and stitching: converting sets of polygons to manifold surface[END_REF]) or perturbations ( [START_REF] Rossignac | Matchmaker: manifold brep for nonmanifold r-sets[END_REF]), but the result is not robust to bad points since the geometry is (almost) the same. Thus previous sparse methods ( [START_REF] Faugeras | Representing stereo data with the delaunay triangulation[END_REF][START_REF] Pan | Proforma: probabilistic feature-based online rapid model acquisition[END_REF][START_REF] Jancosek | Multi-view reconstruction preserving weakly-supported surfaces[END_REF][START_REF] Vu | High accuracy and visibility-consistent dense multiview stereo[END_REF][START_REF] Hoppe | Incremental surface extraction from sparse structure-from-motion point clouds[END_REF][START_REF] Sugiura | 3d surface extraction using incremental tetrahedra carving[END_REF][START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF][START_REF] Litvinov | Incremental solid modeling from sparse structure-from-motion data with improved visual artifact removal[END_REF][START_REF] Piazza | Real-time cpu-based large-scale three-dimensional mesh reconstruction[END_REF]) estimate another bipartition of T , that is the first one up to a regularization. Sec. 1.3.1 will summarize these works if they enforce the manifold constraint.

Other sparse methods ( [START_REF] Morris | Image-consistent surface triangulation[END_REF][START_REF] Vogiatzis | Bayesian stochastic mesh optimization for 3d reconstruction[END_REF][START_REF] Taylor | Surface reconstruction from feature based-stereo[END_REF][START_REF] Hilton | Scene modeling from sparse 3d data[END_REF][START_REF] Lhuillier | A generic error model and its application to automatic 3d modeling of scenes using a catadioptric camera[END_REF][START_REF] Bodis-Szomoru | Efficient edge-aware surface mesh reconstruction for urban scenes[END_REF]) use 2D (not 3D) Delaunay triangulations constrained in images by interest points and image contours. However, [START_REF] Morris | Image-consistent surface triangulation[END_REF][START_REF] Vogiatzis | Bayesian stochastic mesh optimization for 3d reconstruction[END_REF] are limited to simple topology (plane/sphere), large scale scenes cannot be computed efficiently by [START_REF] Taylor | Surface reconstruction from feature based-stereo[END_REF] due to the use of regular decomposition of the 3D space, and [START_REF] Hilton | Scene modeling from sparse 3d data[END_REF][START_REF] Lhuillier | A generic error model and its application to automatic 3d modeling of scenes using a catadioptric camera[END_REF][START_REF] Bodis-Szomoru | Efficient edge-aware surface mesh reconstruction for urban scenes[END_REF] do not generate a globally consistent manifold surface but several (view-centered) local models with redundancy. There are also methods that need strong assumptions to segment the scene using planar regions ( [START_REF] Werner | New techniques for automated architecture reconstruction from photographs[END_REF][START_REF] Fraundorfer | Piecewise planar scene reconstruction from sparse correspondences[END_REF]) or swept surfaces for architectural scenes ([43]).

Previous works enforcing topology constraints

Sec. 1.3 excludes works that are far from our context, e.g. the dense stereo refinements that deform a surface without topology change. (Number of connected components and genus are fixed.)

Manifoldness

Here we focus on the surface reconstruction methods with same input as ours: a sparse point cloud estimated from images and their visibility information. Works ( [START_REF] Faugeras | Representing stereo data with the delaunay triangulation[END_REF][START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF][START_REF] Litvinov | Incremental solid modeling from sparse structure-from-motion data with improved visual artifact removal[END_REF][START_REF] Piazza | Real-time cpu-based large-scale three-dimensional mesh reconstruction[END_REF]) first compute a 3D Delaunay triangulation T of the points and the free-space/matter bipartition of T (Sec. 1.2), then use the manifold property as a constraint to search a surface interpolating most points. In works ([22, 25, 32]), a second bipartition of T in outside tetrahedra and their inside complement is computed by growing an outside set O in the free-space: O should be as large as possible such that the set of the triangles separating outside and inside (i.e. the boundary of O) forms a manifold surface. Visual artifacts of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] are corrected by [START_REF] Litvinov | Incremental solid modeling from sparse structure-from-motion data with improved visual artifact removal[END_REF] and are limited by [START_REF] Piazza | Real-time cpu-based large-scale three-dimensional mesh reconstruction[END_REF], but the artifact problem is not completely solved ( [START_REF] Lhuillier | Overview of shelling for 2-manifold surface reconstruction based on 3d delaunay triangulation[END_REF]). The visibility is partly ignored by [START_REF] Faugeras | Representing stereo data with the delaunay triangulation[END_REF] since the inside grows in the matter.

We note that methods ( [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF][START_REF] Litvinov | Incremental solid modeling from sparse structure-from-motion data with improved visual artifact removal[END_REF][START_REF] Piazza | Real-time cpu-based large-scale three-dimensional mesh reconstruction[END_REF]) are shrink-wrapping methods in the discrete space T guided by the visibility and the manifold constraint. Previous shrink-wrapping methods ( [START_REF] Kobbelt | A shrink wrapping approach to remeshing polygonal surfaces[END_REF][START_REF] Hagbi | Carving for topology simplification of polygonal meshes[END_REF]) ignore the visibility, [START_REF] Hagbi | Carving for topology simplification of polygonal meshes[END_REF] do not use the manifoldness as a constraint on the surface since they use the method of [START_REF] Gueziec | Cutting and stitching: converting sets of polygons to manifold surface[END_REF], the surface has zero-genus by [START_REF] Kobbelt | A shrink wrapping approach to remeshing polygonal surfaces[END_REF].

Low genus

The topological noise is usually removed, or the genus decreases, by a post-processing of the surface reconstruction: [START_REF] Wood | Removing excess topology from isosurfaces[END_REF][START_REF] Zhou | Topology repair of solid models using skeletons[END_REF] remove spurious handles (assuming that they are smaller than true handles) and [START_REF] Hagbi | Carving for topology simplification of polygonal meshes[END_REF] provide multiresolution surfaces including topology simplifications. Now we focus on methods enforcing a low genus during the surface reconstruction. This has been reported in very few literature.

Assuming that the surface topology is known (both number of connected components and genus), [START_REF] Zou | Topology-constrained surface reconstruction from cross-sections[END_REF] show that it is possible to reconstruct an accurate surface from a few planar crosssections. Note that a cross-section is nothing but a sparse cloud of points forming closed curves in R 3 .

In a different context, [START_REF] Hornung | Robust reconstruction of watertight 3d models from non-uniformly sampled point clouds without normal estimation[END_REF] estimate a surface by using a regular grid of voxels (which is inefficient for large scale scenes like ours) to discretize the unsigned distance function to the input points. Compared to the other volumetric methods that use a signed distance function, this method reduces topological noise artifacts caused by misalignments of 3d scan points.

Our previous work ( [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF]) removes the largest ("visually critical") handles generated by a step of the method. Such a handle removal contributes to the visibility consistency optimization, since the outside O increases in the free-space. (It is not like a post-processing ignoring the visibility.) At first glance this quantitatively enforces low genus during surface reconstruction, but we show that this is not the case.

Connectedness

The connectedness is enforced by [START_REF] Oswald | Generalized connectivity constraints for spatio-temporal 3d reconstruction[END_REF] to preserve fine scale details (such as a stick or a rope) without decreasing the genus. First an approximate visual hull including the surface is computed thanks to a green room. Then a convex optimization problem is solved: find an occupancy labeling (sampled in a regular grid) that minimizes a photo-consistency score subject to linear constraints on labeling derivative. The constraints are defined from the visual hull such that the final surface and the visual hull boundary have similar topology.

In contrast to the previous example which enforces inside connectedness, we use outside connectedness like [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF][START_REF] Litvinov | Incremental solid modeling from sparse structure-from-motion data with improved visual artifact removal[END_REF]: the outside O grows continuously in the free-space during most computation. This constraint is legitimate assuming continuous camera motion and opaque matter in the scene.

Our contributions

The surface reconstruction method of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] is the basis of our work. We improve it in different ways: acceleration thanks to a new manifold test (Sec. 3), more efficient repair for singularity removal (Sec. 4), enforce low genus (Sec. 5), more efficient escapes from local extrema (Sec. 6). Only a part of Sec. 5 was published by [START_REF] Litvinov | Incremental solid modeling from sparse structure-from-motion data with improved visual artifact removal[END_REF]. Now we describe the four contributions.

The new manifold test is important since it is used in almost all operations of our methods. It is faster than the previous test by [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] for reasons that were ignored before: the 3D Delaunay triangulation is orientable and the number of surface triangles including a vertex is small in practice. We also provide details on the orientability itself since this property (in the 3D case) is rarely described in the surface reconstruction bibliography.

The new repair operation succeeds more often than the old repair by [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] since the former is guided by an analysis of the surface singularities. (The latter is not.) The former is also quite faster; the latter is the most time consuming operation used by [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF].

Euler's formula (unused by [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF][START_REF] Litvinov | Incremental solid modeling from sparse structure-from-motion data with improved visual artifact removal[END_REF]) reveals that the genus of the O boundary increases a lot due to operations of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF]. We present two methods that greatly reduce the genus. The main idea is to apply such an operation only if it generates a visually non-negligible update of the O boundary at the camera locations used to reconstruct the input points. Up to our knowledge, these methods are the only ones that simultaneously enforce low genus and visibility consistency in surface reconstruction.

Local extrema are partly due to limitations of an operation called "Shelling", which adds one tetrahedron at once to the growing set O. Although it provides most tetrahedra in O, it can "get stuck" in unexpected cases ( [START_REF] Lhuillier | Overview of shelling for 2-manifold surface reconstruction based on 3d delaunay triangulation[END_REF][START_REF] Ziegler | Shelling polyhedral 3-balls and 4-polytopes[END_REF]). We propose a method based on the new repair to escape from local extrema.

In the experiments (Sec. 7), we check the four successive improvements provided by Secs. 3, 4, 5 and 6. We also compare our final method to a graph-cut method by [START_REF] Vu | High accuracy and visibility-consistent dense multiview stereo[END_REF]. Here we do not focus on incremental surface reconstruction like [START_REF] Litvinov | Incremental solid modeling from sparse structure-from-motion data with improved visual artifact removal[END_REF][START_REF] Piazza | Real-time cpu-based large-scale three-dimensional mesh reconstruction[END_REF], although our contributions can be applied to it. (Incremental methods have their own topics.)

Prerequisites

Main notations and definitions

Our input is a set of 3D points and their visibility, i.e. every point p i ∈ R 3 is reconstructed from viewpoints c j ∈ R 3 such that j ∈ V i . A line segment p i c j is called a ray if j ∈ V i . Let T be the 3D Delaunay triangulation computed from the p i ; T is a tetrahedron set that discretizes the convex hull of the p i . The rays define a bipartition of T : the set F ⊆ T of the free-space tetrahedra that are intersected by at least one ray, and the others (the matter tetrahedra).

A simplex σ is the convex hull of k

+ 1 points v 0 • • • v k in general position in R 3 , i.e. v 1 -v 0 • • • v k -v 0 are linearly inde- pendent.
If k is 0, 1, 2 or 3, σ is a vertex, edge, triangle, or tetrahedron, respectively. A simplex σ is a face of σ if σ is the convex hull of some of the v i above. (Thus σ ⊆ σ.) Two tetrahedra are adjacent if they have a common triangle face. Vertices have bold fonts, e.g. v is a vertex and ab is an edge. We use the following notations: ∆ is always a tetrahedron, the sets of the tetrahedra in T that include v or ab are

T v = {∆ ∈ T, v ∈ ∆} and T ab = {∆ ∈ T, ab ⊆ ∆}. ( 1
)
If W is a set of simplices, c(W) is the closure of W, i.e. the set of the faces (including the vertices) of the simplices in W. We have W ⊆ c(W). If K ⊆ c(T ), we sometimes write K instead of the union |K| of its simplices, e.g. we say that K is connected.

If X ⊆ T , ∂X is the boundary of X, i.e. the set of every triangle that is a face of exactly one tetrahedron in X. We have ∂X ⊆ c(X) and every triangle in c(X) \ ∂X is a face of exactly two tetrahedra in X (one in each side of the triangle).

Let v ∈ c(∂X), i.e. v is a vertex of a triangle in ∂X. We say that v is regular in ∂X if the triangles in ∂X including v form a ring around v: the set {ab, abv ∈ ∂X} is a cycle ( [START_REF] Boissonnat | Geometric structures for three-dimensional shape representation[END_REF][START_REF] Lhuillier | 2-manifold tests for 3d delaunay triangulation-based surface reconstruction[END_REF]). Otherwise, v is singular in ∂X. Furthermore, ∂X is manifold iff (if and only if) every vertex in c(∂X) is regular in ∂X. In general, ∂F is not manifold. Let O be a set of outside tetrahedra: we have O ⊆ F ⊆ T and ∂O is manifold. Let g be the genus of ∂O.

Operations on the set of outside tetrahedra

Here we summarize operations introduced by [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF][START_REF] Litvinov | Incremental solid modeling from sparse structure-from-motion data with improved visual artifact removal[END_REF]. They are basic components of methods considered in this paper. Technical details are omitted (as most as possible) to make easier the paper understanding.

Principles

The target of our methods is O ⊆ F that maximizes a visibility score function f (O) subject to the constraint that ∂O is manifold. For efficiency, f (O) is defined as the sum for each tetrahedron ∆ ∈ O of a positive real f (∆). We choose f (∆) as the number of ray intersection(s) for ∆, but alternative definitions of f could be investigated inspired by [START_REF] Pan | Proforma: probabilistic feature-based online rapid model acquisition[END_REF][START_REF] Piazza | Real-time cpu-based large-scale three-dimensional mesh reconstruction[END_REF] ) From top to bottom: Shelling (Sh.), Topology Extension (T.E.), Force-Repair (F.R.) and Shrink-Grow (S.G.). Every line shows successive values of O in gray. Sh. adds one tetrahedron at once to O and g does not change. T.E. adds all tetrahedra around a vertex to O so that g can increase. F.R. adds tetrahedra, a singularity is first created and then removed. S.G. first removes then adds tetrahedra to O. be seen as cross-sections of the former, e.g. the cross-section inside a torus is an annulus or a disconnected set of triangles.) Shelling (Sh.) adds one tetrahedron at once to O such that ∂O is always manifold and the number of connected component(s) of O does not change. The first line of Fig. 2 shows an example. In more details, Sh. builds a tetrahedron series

O 0 , O 1 , • • • O n such that O i+1 = O i ∪ {∆ i } and the tetrahedron ∆ i ∈ F \ O i has at least one triangle face in ∂O i (if O i ∅)
and every ∂O i is manifold. We first consider ∆ i candidates that have the largest f (∆ i ), such that O globally grows from the most confident freespace to the less confident one. Once ∆ 0 is selected, the other ∆ i are selected in a growing neighborhood of ∆ 0 thanks to a priority queue. In practice, Sh. is fast and adds almost all tetrahedra in O. A local Shelling is possible by specifying the first added tetrahedron ∆ 0 as input. More details can be found in algorithm 1 of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF].

However, Sh. is not sufficient since it cannot change g ( [START_REF] Lhuillier | Overview of shelling for 2-manifold surface reconstruction based on 3d delaunay triangulation[END_REF]), the genus of ∂O. If the genus of the true surface is not 0, e.g. if the camera path (the c j series) includes a loop around a high building, g should increase.

Topology Extension (T.E.) builds

O 1 = O 0 ∪ T v if a vertex v ∈ c(∂O 0 ) and T v ⊆ F. It is successful if ∂O 1 is
manifold. Now g can increase as shown in the second line of Fig. 2,since T.E. adds to O several tetrahedra at once. However, T.E. can also generate spurious handles as shown in the left of Fig. 15. Thus another operation is needed; it takes as input a tetrahedron set H ⊆ F \O that we would like to add to O. For example, H is a cross section of a handle that we would like to remove.

Force-and-Repair (F.R.) builds a series O 0 , O 1 , • • • O n defined as follows. First we "force" O 1 = O 0 ∪ H and ignore singular vertices in ∂O 1 . Then we try to "repair". Let s i be the number of singular vertices of ∂O i . We select a tetrahedron ∆ i ∈ F \ O i having triangle face(s) in ∂O i and set O i+1 = O i ∪ {∆ i } while s i+1 ≤ s i . F.R. is successful if s n = 0. More details on "repair" can be found in Algorithm 3 of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] (using notation G = H).

The third line of Fig. 2 shows a F.R. example.

Steps Critical Handle Removals (C.H.R.s) of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] and [START_REF] Litvinov | Incremental solid modeling from sparse structure-from-motion data with improved visual artifact removal[END_REF] select several sets H and try to add each of them to O by using F.R. The C.H.R.s remove spurious handles (Sec. 1.3.2) that are visually important. More details can be found in Algorithm 2 of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] and Sec. III.C of [START_REF] Litvinov | Incremental solid modeling from sparse structure-from-motion data with improved visual artifact removal[END_REF].

Now we note that all operations above (Sh., T.E., F.R.) and their compositions (C.H.R.s) are O growing in the free-space F. By combining them, we obtain a descent method for minimizing functionf (O) such that ∂O is manifold, which can get stuck to a local extremum. The following operation removes some tetrahedra from O (Thus f (O) decreases temporarily.) to kick the algorithm out of local extrema.

Shrink-and-Grow (S.G.) builds a series that is not 100% growing as shown in the last line of Fig. 2. First O shrinks:

O 1 = O 0 \ T v if a vertex v ∈ c(∂O 0 ) and ∂O 1 is manifold. Then O grows: we try local Shelling O 1 , O 2 , • • • O n started by a tetra- hedron ∆ ∈ F \ O 0 that has the vertex v. S.G. is successful if f (O n ) ≥ f (O 0 ).
In practice, S.G. is tried at several selected vertices v ∈ c(∂O). More details can be found in algorithm of Sec. III.B of [START_REF] Litvinov | Incremental solid modeling from sparse structure-from-motion data with improved visual artifact removal[END_REF]. (See also Appendix D.) In the paper remainder, we only need this composition of several S.G. operations and directly use notation S.G. for the composition.

Directed edge-based manifold test

Manifold test at a vertex is important since it is used in all operations in Sec. 2.2.2 except Shelling. (Shelling has its own manifold test.) Sec. 3.1 and 3.2 are reminders about simplex orientation and orientability of T . Then Sec. 3.3 introduces a method that checks if a vertex is regular in ∂O. It is based on a theorem, whose proof is in Sec. 3.4.

Simplex orientation

According to [START_REF] Zomorodian | Topology for computing[END_REF], an orientation of a simplex

v 0 v 1 • • • v k is an equivalence class of its vertex orderings such that (v 0 , v 1 , • • • v k ) and (v π0 , v π1 , • • • v πk )
are equivalent iff the permutation π is even (using shortened notation πi = π(i)). We remind that every permutation is a product of transposition(s), and a permutation is even if it is the product of an even number of transpositions. Furthermore, if there are distinct integers a and b such that πa = b and πb = a and πc = c where c ∈ N \ {a, b}, π is a transposition and we use notation π = (a b).

Let A k+1 be the set of the even permutations of

{0, • • • k}. Thus (v 0 , v 1 , • • • , v k ) = {(v π0 , v π1 , • • • v πk ), π ∈ A k+1 } (2) is the equivalence class of (v 0 , v 1 , • • • v k ) and we have π ∈ A k+1 ⇐⇒ (v 0 , v 1 , • • • v k ) = (v π0 , v π1 , • • • v πk ). (3)
Here are examples using π 1 = (0 1)•(2 3) and π 2 = (0 2)•(0 1): We note that there are exactly two possible orientations if

(0 1) A 4 ⇒ (v 0 , v 1 , v 2 , v 3 ) (v 1 , v 0 , v 2 , v 3 ) (4) π 1 ∈ A 4 ⇒ (v 0 , v 1 , v 2 , v 3 ) = (v 1 , v 0 , v 3 , v 2 ) (5) π 2 ∈ A 4 ⇒ (v 0 , v 1 , v 2 , v 3 ) = (v 1 , v 2 , v 0 , v 3 ) (6) (v 0 , v 1 , v 2 ) = (v 1 , v 2 , v 0 ) = (v 2 , v 0 , v 1 ). (7)
k ≥ 1: (v 0 , v 1 , v 2 , • • • v k ) and (v 1 , v 0 , v 2 , • • • v k ).
Here we swap v 0 and v 1 as in Eq. 4, but other vertices can be swapped. We also note that this orientation definition does not require vertex coordinates, although it is related to orientation in vector

spaces: function (v 0 , • • • , v k ) → det v 0 • • • v k 1 • • • 1 is constant
in an equivalence class and changes its sign in the other class.

Orientability of T

We first provide the definition of induced orientation since it helps to understand orientability. An orientation

(v 0 , v 1 , • • • v k ) induces orientations of faces of simplex v 0 v 1 • • • v k as follows: the induced orientation of v π1 • • • v πk is (v π1 , • • • v πk ) if π ∈ A k+1 . For example, a triangle orientation (v 0 , v 1 , v 2 ) induces edge ori- entations (v 1 , v 2 ), (v 2 , v 0 ) and (v 0 , v 1 ) thanks to Eq. 7. Similarly, a tetrahedron orientation (v 0 , v 1 , v 2 , v 3 ) induces triangle orien- tations (v 1 , v 2 , v 3 ), (v 2 , v 0 , v 3 ), (v 0 , v 1 , v 3 ) and (v 2 , v 1 , v 0 )
. The left of Fig. 3 shows these induced orientations by using a standard convention: the orientation (v a , v b , v c ) is represented by an arrow that rotates in the same direction as

v a → v b → v c → v a .
Now assume that we have orientations (v 0 , v 1 , v 2 , v 3 ) and (v 0 , v 1 , v 2 , v 3 ) of adjacent tetrahedra ∆ and ∆ , respectively. We say that ∆ and ∆ are consistently oriented if there are π ∈ A 4 and π [START_REF] Teillaud | Three dimensional triangulation in cgal[END_REF]). Fig. 3 shows tetrahedra abcd and abce with consistent orienta- tions (a, b, c, d) and(b, a, c, e). We check1 that they induce different orientations (b, a, c) and (a, b, c) for their shared triangle face abc. Intuitively, d and e "see" abc in opposite orientations.

∈ A 4 such that (v π0 , v π1 , v π2 ) = (v π 1 , v π 0 , v π 2 ) ([
Last T is orientable ( [START_REF] Teillaud | Three dimensional triangulation in cgal[END_REF][START_REF] Boissonnat | Triangulations in cgal[END_REF], see also Proofs.pdf in the supplementary material), i.e. Theorem 1. There are orientation choices for all tetrahedra in T such that every pair of adjacent tetrahedra in T is consistently oriented.

Manifold test

We first introduce notations and definitions to check that a vertex v ∈ c(∂O) is regular in ∂O. Theorem 1 provides an orientation o(∆) for every tetrahedron ∆ ∈ T . This is done once before all manifold tests. We consider a set of directed edges:

D = ∆∈O {(v 2 , v 3 ), (v, v 1 , v 2 , v 3 ) ∈ o(∆), vv 2 v 3 ∈ ∂O}. (8)
If D = {(q 1 , q 2 ), • • • (q m-1 , q m ), (q m , q 1 )} and m ≥ 3 and

q 1 • • • q m are m distinct vertices, D is a directed cycle.

Theorem

Our test is based on the following theorem.

Theorem 2. The vertex v ∈ c(∂O) is regular in ∂O iff D is a directed cycle.
We check Theorem 2 in simple examples to provide intuition. In the first example, a tetrahedron vabc ∈ O and o(vabc) = (v, a, b, c) and triangles vbc, vca, vab are in ∂O. Since

{(v, a, b, c), (v, b, c, a), (v, c, a, b)} ⊂ o(vabc), (9)

Implementation

Here we summarize the implementation of the manifold test. For every vertex w ∈ c(T ), T w is stored in a table. We store D in a table E of vertices where the i-th directed edge of D has its start-vertex in E[2i] and its end-vertex in E

[2i+1]. For every tetrahedron ∆ = v ∆ 0 v ∆ 1 v ∆ 2 v ∆ 3 ∈ T , we store a vertex ordering (v ∆ 0 , v ∆ 1 , v ∆ 2 , v ∆ 3 ) such that o(∆) = (v ∆ 0 , v ∆ 1 , v ∆ 2 , v ∆ 3 ). Thus, D = ∆∈O {(v ∆ π2 , v ∆ π3 ), π ∈ A 4 , v ∆ π0 = v, v ∆ π0 v ∆ π2 v ∆ π3 ∈ ∂O}. ( 10
)
Thanks to Eq. 10, for every

∆ ∈ T v ∩ O and π ∈ A 4 , we collect in E the edges using a C-like line as if (v ∆ π0 ==v && v ∆ π0 v ∆ π2 v ∆ π3 ∈ ∂O) { E[n++]=v ∆ π2 ; E[n++]=v ∆ π3 ;
} after initializing n=0. Last we check that D is a directed cycle using a naive sorting algorithm: swap directed edges in the table [START_REF] Boissonnat | Triangulations in cgal[END_REF] and so on.

E to form path E[0]E[1], E[2]E[3] with E[1]==E[2], E[4]E[5] with E[3]==E
Appendix A gives more details about the algorithm. This appendix also details a previous test algorithm of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] that is compared in our experiments. The efficiency of the new test is based on facts that are not considered before: small D size and consistent orientations in T . Thus the naive sorting in E is fast in spite of its complexity (quadratic in the size of D).

Proof of Theorem 2

This paragraph only provides the principle of the proof; the details are in Appendix B. The proof may look difficult at first glance, but it should be stressed that the corresponding algorithm is not: Algorithm 3 only reads directed edges and naively sorts them in a table to detect a directed cycle.

The undirected version of D is

U = {ab, (a, b) ∈ D}. (11)
We note that U is a cycle (with undirected edges) if D is a directed cycle (with directed edges), but the converse can be false. First it is not difficult to see that Lemma 1. U is the set of the v-opposite edges in the triangles of ∂O:

U = {ab, abv ∈ ∂O}. Since v is regular iff {ab, abv ∈ ∂O} is a cycle, Lemma 2. The vertex v is regular in ∂O iff U is a cycle.
The easy part of Theorem 2's proof is the sufficient condition "D is a directed cycle ⇒ v is regular" thanks to Lemma 2, which simplifies this to "D is a directed cycle ⇒ U is a cycle", which in turn is straightforward. The necessary condition "v is regular ⇒ D is a directed cycle" is more difficult although Lemma 2 simplifies this to "U is a cycle ⇒ D is a directed cycle". The main part of the proof is to show that the consistent orientations in T imply good edge directions in D; this is done by studying consistent orientations of the tetrahedra sharing a common edge face vu ∈ c(∂O).

Repair by analyzing the surface singularities

Sec. 4.3 presents a new repair operation that is guided by an analysis of the surface singularities to be removed. The analysis is done by using two Theorems in Sec. 4.2, which count the connected components of adjacency graphs of tetrahedra. We first extend T (Sec. 4.1) to avoid special cases in the Theorems.

Abstract extension T ∞ of T

To avoid special cases in statements and proofs, i.e. if c(∂T ) includes singular vertices of ∂O, we extend T into T ∞ such that ∂T ∞ = ∅ like [START_REF] Boissonnat | Triangulations in cgal[END_REF][START_REF] Lhuillier | 2-manifold tests for 3d delaunay triangulation-based surface reconstruction[END_REF].

First we replace every vertex of c(T ) by an integer that identifies the vertex. Then every simplex of c(T ) with 2,3,4 vertices is replaced by a set of 2,3,4 integers that identifies its vertices, respectively. (e.g. an edge becomes the set of the two integers of its vertices.) The relation "is a face of" (inclusion) does not change. Bold fonts and words vertex/edge/triangle/tetrahedron are also used for this new "abstract" simplex definition.

Then we extend T . Let v ∞ be an integer that is different to all vertices of c(T ). For every triangle abc ∈ ∂T , we create a new tetrahedron abcv ∞ by adding v ∞ to integer set abc. Let

T ∞ = T ∪ {abcv ∞ , abc ∈ ∂T }. ( 12
)
Now T ⊂ T ∞ and every triangle in c(T ∞ ) is a face of exactly two tetrahedra in T ∞ . Fig. 4 shows v ∞ and T ∞ . The tetrahedron sets in T ∞ that include a vertex v or an edge ab in c(T ∞ ) are

T ∞ v = {∆ ∈ T ∞ , v ∈ ∆} and T ∞ ab = {∆ ∈ T ∞ , ab ⊆ ∆}. ( 13
)

Vertex singularities and edge singularities

Let vertex v ∈ c(∂O) and g v be the adjacency graph of the tetrahedra in T ∞ v . (There is a bijection between the vertices of g v and the tetrahedra in T ∞ v , every edge of g v links two vertices if their corresponding tetrahedra are adjacent.) Let g O v be the graph obtained from g v by removing every edge between a tetrahedron in T ∞ v ∩ O and a tetrahedron in

T ∞ v \ O. Thus Lemma 3. Every connected component of g O v is included in T ∞ v ∩ O or included in T ∞ v \ O. There is at least one component in T ∞ v ∩ O and at least one component in T ∞ v \ O. Proof. Since all edges between T ∞ v ∩O and T ∞ v \O are removed from g O v , the first assertion is true. There is a triangle t ∈ ∂O including v; t is a face of a tetrahedron in T ∞ v ∩ O and is a face of another one in T ∞ v \ O.
Thus the second assertion is true.

Theorem 3. The vertex v is singular in ∂O iff g O v has at least three connected components. (See Fig. 5, see also Fig. 4 of [START_REF] Lhuillier | 2-manifold tests for 3d delaunay triangulation-based surface reconstruction[END_REF].)

Proof. Thanks to Lemma 3, g O v has at least two connected components. Thanks to Theorem 3 of [START_REF] Lhuillier | 2-manifold tests for 3d delaunay triangulation-based surface reconstruction[END_REF], v ∈ c(∂O) is singular iff g O v has at least three connected components.

Similarly, let edge ab ∈ c(∂O) and g ab be the adjacency graph of the tetrahedra in T ∞ ab . Let g O ab be the graph obtained from g ab by removing every edge between a tetrahedron in T ∞ ab ∩ O and a tetrahedron in T ∞ ab \ O. Thus Lemma 4. The graph g ab is a cycle. The graph g O ab has exactly 2n connected components: n in T ∞ ab ∩O and n in T ∞ ab \O (Fig. 6). Proof. Since T ∞ is an (extended) 3D Delaunay triangulation,

T ∞ ab = {abc 0 c 1 , abc 1 c 2 , abc 2 c 3 , • • • abc k c 0 } ( 14
)
where the edges

c 0 c 1 , • • • c k-1 c k , c k c 0 form a cycle. Thus g ab is a cycle. Furthermore, all edges between T ∞ ab ∩ O and T ∞ ab \ O are removed in g O ab . Thus every connected component of g O ab is a path in T ∞ ab ∩ O or in T ∞ ab \ O. Since every path in T ∞ ab ∩ O is before a path in T ∞
ab \O in the cycle (or conversely), the numbers of components in T ∞ ab ∩ O and T ∞ ab \ O are the same.

We say that the edge ab is singular in ∂O if g O ab has at least three connected components. Thus Theorem 4. Assume that the edge ab is singular in ∂O. Then g O ab has at least four connected components (Fig. 6). Furthermore, both vertices a and b are singular in ∂O.

Proof. Thanks to Lemma 4, the first assertion is true. For the second assertion, we assume that a is regular in ∂O and show that g O ab has exactly two connected components. Thanks to Lemma 4, ∂O has exactly 2n triangles including ab. (Such a triangle is between T ∞ ab ∩ O and T ∞ ab \ O.) Since a is regular, there are exactly two triangles in ∂O having the edge ab ( [START_REF] Lhuillier | 2-manifold tests for 3d delaunay triangulation-based surface reconstruction[END_REF]). Thus n = 1 and g O ab has exactly two connected components.

Method

Singular vertices appears while adding to O a tetrahedron set H ⊆ F \ O. Our repair tries to remove them by adding to O other tetrahedra in F \ O without creating new singular vertices.

Principle

If ab is a singular edge of ∂O, g O ab has at least two connected components in T ∞ ab \ O (right of Fig. 6). We can add to O such a connected component if it is included in F and no additional singular vertex appears. Then the number of connected components of g O ab decreases (middle of Fig. 6), and we progress toward a non-singular ab thanks to Theorem 4. If ab becomes non-singular, vertices a or/and b can become non-singular too.

If a is a singular vertex of ∂O, g O a has at least one connected component in T ∞ a \ O (e.g. two on the right of Fig. 5). Similarly, we can add to O such a connected component if it is included in F and no additional singular vertex appears. Then the number of connected components of g O a decreases (left of Fig. 5) and we progress toward a non-singular a thanks to Theorem 3. If T ∞ a \O becomes empty, a become non-singular since a c(∂O).

At first glance, the adds near singular edges are useless since

• a singular edge is a special case of a pair of singular vertices according to Theorem 4 and

• ∂O is manifold iff its vertices are non-singular (Sec. 2.1). This is wrong: there are cases where a singular edge ab is only removable by an add near the edge since the adds near vertices a and b are impossible. For example, assume that T ∞ ab has four tetrahedra:

abc 0 c 1 ∈ T \ F, abc 1 c 2 ∈ O,
i = true && b 1 i = f alse) { O ← O \ C; return 0; } 07: for each i ∈ {0, 1, • • • c} 08: if (b 0 i = f alse && b 1 i = true) V ← V \ {v i }; 09: A ← A ∪ C; 10: return 1; // success end
Then we detail our Repair2 method in Algorithm 2 (in C style). It takes as input O and the set H that was just added to O. (∂O was manifold before that.) It use a set E where the singular edges are searched. It also use a set A that collects the new tetrahedra added to O; A is used to restore O in case of failure or to improve O in case of success. Repair2 examines singular edges (special cases) before singular vertices. The set V can only decrease and Repair2 is successful iff V = ∅. Fig. 7 shows an example of Repair2. The main loop is not infinite, since its number of iterations is less than the number of tetrahedra in F \ O having a vertex in c(H). In our implementation, the sets H, V, E, C and A are stored in tables. end Last we explain why Repair (Algorithm 3 of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF]) is less efficient than Repair2 above. Repair ignores the structure of the singularities (in terms of connected components of adjacent tetrahedra) and does not distinguish singular edges from general singular vertices; it blindly adds tetrahedra one-by-one to O such that no new singular vertices appears. Furthermore, Repair is time consuming since it applies manifold tests before and after adding every tetrahedron. Thanks to V, Repair2 only needs to apply manifold tests after adding packs of tetrahedra.

Algorithm 2. Repair2

E 11 if (a ∈ V && b ∈ V) { 12: Let C 1 , • • • , C k be the connected components of g O ab ; 13: if (there is a C i such that C i ⊆ F \ O && 14: ReduceSingularity(C i , V, O, A)) IsImproved= 1; 15: } else E ← E \ {ab}; 16: for each vertex a ∈ V { 17: Let C 1 , • • • , C k be the connected components of g O a ; 18: if (there is a C i such that C i ⊆ F \ O

Lowered genus during surface reconstruction

Sec. 5.1 shortly summarizes the method of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] and explains how the genus g of ∂O evolves during the computation. Sec. 5.2 and 5.3 present two methods estimating ∂O with a smaller g. We remember Euler's formula ( [START_REF] Botsch | Polygon mesh processing[END_REF]) if ∂O is connected:

2(1 -g) = v -e + t (15)
where v, e and t are the numbers of vertices, edges and triangles in c(∂O), respectively. If ∂O is not connected, its genus is the sum of genuses of its connected components.

Previous method and g evolution

The previous method of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] has several steps defined from operations summarized in Sec. 2.2.2.

Step 1: Shelling plus Topology Extension

The set O is initialized by Sh., which tries to add first to O the tetrahedra that have the largest number of ray intersections. (The first tetrahedron in O maximizes f .) The resulting ∂O meets g = 0 although the true surface has a non-zero genus, e.g. if the camera path has loops around buildings. Then T.E. is tried on every vertex v ∈ c(∂O). Now g > 0, but spurious handles occur.

Step 2: Critical Handle Removal

C.H.R. in Sec. 4.3 of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] is applied to remove critical handles, i.e. spurious handles that the human eye cannot miss at the camera locations c j where the images are taken. Here we provide useful details on C.H.R.

First we define critical edges ab for a given angle threshold α > 0: ab is a face of a tetrahedron in F \ O (which can be added to O), it is only included in tetrahedra in F, it is large enough such that angle ac j b > α. Thus the set of the critical edges is concisely defined by

L α = {ab ∈ c(F \ O) \ c(∂T ), T ab ⊆ F, ∃ j, ac j b > α}. ( 16
)
Then C.H.R. subdivides every critical edge by adding a Steiner vertex at the middle to create new tetrahedra, such that ∂O is still manifold. (A tetrahedron can be split in two new tetrahedra that are in O or F iff the original one is.) A Steiner vertex is an extra point that is not in the input cloud. Last F.R. tries to add to O a tetrahedron set H included in T v ∩ (F \ O), for every vertex v of the split edges.

Although critical handles are removed like this, our experiment shows that g increases a lot during this step.

Step 3: post-processing including Peak Removal

The post-processing improves ∂O thanks to prior knowledge of the scene. It includes a step Peak Removal (P.R), a Laplacian smoothing and a process to remove triangles in the sky. The two latters are not summarized since they do not change g. P.R. removes bad points in the surface. It is tried on every vertex v ∈ c(∂O) such that the ring of the triangles in ∂O including v has a solid angle w in the T ∞ v \ O side such that w < w 0 . Then it add to O the set T v if ∂O remains a manifold. (The process is similar in the T ∞ v ∩ O side.) It looks like T.E. without condition T v ⊆ F. Our experiment shows that g decreases thanks to P.R.

Method 1: cancel C.H.R. tentatives

The first idea that comes to mind is to cancel C.H.R. if g increases. Let δg be the g variation due to a successful F.R. operation during C.H.R. We also would like an efficient computation of δg without a large traversal of ∂O. Since T is implemented by the adjacency graph of the tetrahedra, large traversals in T should be avoided.

Let c be the number of the connected components of ∂O. By summing Eq. 15 for all connected components, we obtain

δv -δe + δt = 2(δc -δg) (17)
where δv, δe, δt and δc are the variations of v, e, t and c due to a successful F.R. operation. Sec. 5.2.1 presents an efficient computation of δvδe + δt. Then Sec. 5.2.2 presents a method to avoid the computation of δc which is not efficient.

J = (∂O \ ∂(O ∪ H)) ∪ (∂(O ∪ H) \ ∂O). ( 18
)
We note that the triangles in ∂O ∩ ∂(O ∪ H)) do not contribute to δt. (They add 0 to δt.) Similarly, the vertices and edges of these triangles do not contribute to δv and δe. Thus we only need a traversal of J to find all the vertices/edges/triangles that contribute to δvδe + δt. Intuitively, we only need a traversal of H to compute δvδe + δt since the boundary change |J| is included in the volume change | H|. Indeed, Appendix C shows that Lemma 5. We have ∂ H = J.

Thanks to a traversal of H, we access to the triangles in ∂ H and count those of them that are in ∂O. We also count those of them that are in ∂(O ∪ H) and obtain δt. Similarly, we access to the vertices and edges in c(∂ H) and count those of them that are in c(∂O) and those that are in c(∂(O ∪ H)), then we obtain δv and δe. This computation is efficient since H is quite smaller than O and O ∪ H and their boundaries.

Avoid inefficient computation of δc

At first glance, we accept a successful F.R. operation iff δg ≤ 0. Here δg is estimated using Eq. 17 and Sec. 5.2.1 and the computation of δc. Unfortunately, the δc computation is not efficient. Indeed, the top row of Fig. 8 shows that we cannot distinguish cases (δg = -1, δc = 0) and (δg = 0, δc = 1) by a traversal of the small set H; a traversal of large set ∂(O ∪ H) should be done instead.

Then we use condition δcδg ≥ 0 instead of δg ≤ 0 for three reasons. First both conditions are simultaneously true or simultaneously false in all cases of Fig. 8. Second, the computation of δc-δg is efficient thanks to Eq. 17 and Sec. 5.2.1. Last δg < 0 "almost" implies δcδg ≥ 0 according to the following theorem (proof in Sec. 5.2.3). Theorem 5. The set H is connected. If O is connected and ∂O is manifold and H ∩O = ∅, then δc ≥ -1 (in most cases δc ≥ 0). Indeed, the assumptions of Theorem 5 are "almost" meet during C.H.R. of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] since the number of connected components of O (Note: this is not c.) is very low in practice.

Sketch proof of Theorem 5

The proof details are in file Proofs.pdf of the supplementary material. First H is connected since all tetrahedra in H share a same vertex and H is obtained by a continuous growing started from H ( [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF]). Second we see that O is a "solid" ( [START_REF] Sakkalis | Representational validity of boundary representation models[END_REF]): it is like one piece of gruyere with one external boundary and k internal boundaries, each of them is the manifold boundary of a connected component K i of R 3 \ (|O| \ |∂O|). (See middle of Fig. 9.) Third we see that H is included in a single K i . (Since H is connected and H ∩ O = ∅.) Fourth, K j is unchanged by the add of H to O if j i. Last there are at least k connected components of ∂O that are unchanged by the add of H to O, and we conclude.

We note that δc ≥ 0 in most cases according to Fig. 9, i.e. iff H is strictly included in K i . (We have δc = -1 iff H = K i .)

Method 2: replace T.E. and C.H.R. operations

This method has the same start (Sh. defined by Algorithm 1 of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF]) and same end (Sec. 5.1.3) as the previous method in Sec. 5.1. T.E. is replaced by the step in Sec. 5.3.2, and C.H.R. is replaced by the step in Sec. 5.3.3. The genus g moderately increases thanks to these two new steps using ideas in Sec. 5.3.1. We moderate the g increase as follows. First, T.E. is replaced by an operation that is applied only if it provides a visually nonnegligible change of ∂O. In our case, this operation removes a critical edge from c(∂O). Second, the C.H.R. of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] is replaced by a new C.H.R. that is more selective on set H of tetrahedra in F \ O that it tries to add to O. Intuitively, H is the set of the tetrahedra in a handle that can be cut by a plane. Furthermore, c(H) also includes a critical edge.

Critical Edge Removal (C.E.R.)

This step removes critical edges from c(∂O). For every edge ab ∈ L α ∩ c(∂O), F.R. (Sec. 2.2.2) is applied using H = T ab ∩ (F \ O). In more details, we first apply O ← O ∪ H, then we try to remove the singular vertices in ∂O using Algorithm 3 of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] with input G = H. If the final ∂O is not manifold, the ab removal fails and we restore O to its initial value.

This step acts as T.E. at locations (critical edges) that are selected by α, i.e. g can increase only if this provides a visually important change of ∂O. Fig. 10 shows that g can increase or even decrease by using C.E.R. Note that O is improved in both cases since f (O) increases. Last we improve the result (without g change) thanks to Sh., i.e. we try to start a local Shelling from every tetrahedron ∆ ∈ F \ O.

Critical Handle Removal (C.H.R.)

A critical handle H meets several conditions ( [START_REF] Litvinov | Incremental solid modeling from sparse structure-from-motion data with improved visual artifact removal[END_REF]). First, H ⊆ F \ O. Second H is "critical": there is a critical edge ab ∈ L α such that ab ∈ c(H). Third, there is a plane π such that π ∩ ab ∅ and π "cuts" H: H ⊆ T π where

T π = {∆ ∈ T, π ∩ ∆ ∅}. Last H is "surrounded" by O in π: we have ∆ ∈ O if ∆ ∈ T π \ H is adjacent to another tetrahedron ∆ ∈ H.
The detection of H is done as follows. For every edge ab ∈ L α , we try several planes π that are orthogonal to ab or horizontal. (These directions are expected for cross-sections of H by π.) Then we search H by growing in T π ∩(F \O) and starting by H = T ab ∩ (F \ O). If the surrounding condition is not meet at the growing end, the detection of H fails and another pair (ab, π) is tried.

Once H is computed, we use F.R. as in C.E.R above: first apply O ← O ∪ H, then try to remove the singular vertices in ∂O thanks to Algorithm 3 of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] using input G = H. Fig. 11 shows an example of C.H.R.

At this point, we can remove remaining critical handles using C.H.R. of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] (summary in Sec. 5.1.2), but we would like to moderate its g increase (due to small handles) and its use of Steiner vertices. Let M be the set of tetrahedra obtained as the union of the remaining handles detected as in the beginning of Sec. 5.3.3. We simply use the C.H.R. summarized in Sec. 5.1.2 by replacing its input set L α by L α ∩ c(M).

Last we note that C.H.R. is preceded by S.G. in [START_REF] Litvinov | Incremental solid modeling from sparse structure-from-motion data with improved visual artifact removal[END_REF]. The summary in Sec. 2.2.2 reminds us that S.G. is useful to escape from local extremum of function f . Since our experiments show that g is almost constant by using S.G., we also use S.G. before C.H.R. (and after C.E.R) in method 2. Details on S.G. are given in Appendix D.

Escape from local extrema due to Shelling blocking

First Sec. 6.1 asserts that Shelling alone cannot generate as many surfaces as expected, then Sec. 6.2 provides examples. Last Sec. 6.3 explains how to solve (partly) this problem.

Definition of Shelling blocking

Let O and O be tetrahedron sets such that andx → h(x, t) is homeomorphism ∀t.) Thus ∂O is manifold with the same genus as ∂O. Let n be the number of the tetrahedra in O \ O. We could expect to obtain O from O by a greedy Shelling algorithm: let

∅ O ⊂ O ⊆ F, ∂O is manifold, |∂O| can be continuously deformed to |∂O | by moving in |O \ O| using an isotopy. (It is a continuous func- tion h : |∂O| × [0, 1] → |O \ O| such that h(|∂O|, 0) = |∂O|, h(|∂O|, 1) = |∂O |,
O 0 = O, then choose tetrahedron series ∆ i ∈ O \ O i-1 for i varying from 1 to n and set O i = {∆ i } ∪ O i-1 (Sec. 2.2.
2). However, this is not always possible. In this case, we say that we have a Shelling blocking.

A Shelling blocking implies that there is i ≤ n such that no tetrahedron ∆ i ∈ O \ O i-1 meets all required constraints to define O i . In other words, ∂(O i-1 ∪ {∆ i }) has a singular vertex or ∆ i has no triangle face in ∂O i-1 for every tried

∆ i ∈ O \ O i-1 .
This may look surprising since Shelling in the 2D case does not have such a blocking according to [START_REF] Danaraj | A representation of 2-dimensional pseudomanifolds and its use in the design of a linear-time shelling processing[END_REF]. (Replace tetrahedra by triangles, 2-manifold by 1-manifold, as in the top line of Fig. 2.)

, O , O \ O is homeomorphic to {x ∈ R 3 , ||x|| ≤ 1}, c(O) ∩ c(O \ O) is a disc.

Examples of Shelling blocking

In a first example, ∂O is a sphere (i.e. ∂O is homeomorphic to {x ∈ R 3 , ||x|| = 1}.) and we set O = {∆} with a tetrahedron ∆ ∈ O . According to [START_REF] Ziegler | Shelling polyhedral 3-balls and 4-polytopes[END_REF][START_REF] Lhuillier | Overview of shelling for 2-manifold surface reconstruction based on 3d delaunay triangulation[END_REF], there are three cases. Case 1: Shelling can success and can fail (This depends on its tetrahedron choices.), e.g. if O is convex and large enough. Case 2: Shelling always fails, whatever its tetrahedron choices. (An example O has only 12 vertices and 25 tetrahedra.) Case 3: there is never Shelling blocking, e.g. if the number of vertices in c(O ) is less than 9.

In a second example, ∂O is a manifold, ∂(O \ O) is a sphere, and c(O) ∩ c(O \ O) is a disc. (i.e. it is homeomorphic to {x ∈ R 2 , ||x|| ≤ 1}.) Fig. 12 shows an example where ∂O is a torus. Appendix E explains how to obtain a Shelling blocking in this context if O \ O has only 5 tetrahedra.

Escape from Shelling blocking using other operations

T.E., C.H.R., F.R., C.E.R. and S.G. in Secs. 2.2.2 and 5.3.2 can escape from Shelling blocking, although they are not originally designed to do that. This means that they can reach O from O if we have a Shelling blocking. (Note that S.G. is not designed to meet O ⊂ O .) Appendix E.1 shows an example: C.E.R. can escape from the last Shelling blocking in Sec. 6.2.

Here we present a method called Unlock to escape from remaining blockings. F.R. is used with two differences compared to [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF]. First Repair is replaced by Repair2 (in Sec. 4.3.2) which is better. Second Repair2 is used without involving visually critical edges to select locations where it is applied. We apply Repair2 to a lot of small H.

(Reminder: first force O ← O ∪ H then repair O in the neighborhood of H.) We try H = {∆} for every tetrahedron ∆ ∈ F \ O and H = T v ∩ (F \ O) for every vertex v ∈ c(F \ O).
However, this generates high topological noise. A solution is to cancel the successful F.R. tentatives as in Sec. 5.2, but this is time consuming.

We prefer to enforce a supplementary constraint on the tried H before F.R.: the graph whose vertices are in c(H) ∩ c(O) and whose edges are in c(∂O) must be non-empty and connected. (The complete edge set is stored in an adjacency list before all F.R. for efficiency.) Experiment shows that this heuristic constraint greatly reduces topological noise and we explain it in an example. Assume that H = {∆} and F \ O is like a coin: a cylinder with large diameter and small thickness (Fig. 13). The constraint avoids almost all cases where ∆ has a vertex in one side and another vertex in the other side of the coin. If F.R. is applied and is successful in such a case (i.e. without the contraint), it creates a hole connecting both sides and the genus of ∂O increases.

Experiments

Sparse input point cloud from images

We reconstruct complete environments thanks to videos taken by helmed-held omnidirectional multi-cameras and spherical cameras.

First the multi-camera is selfcalibrated ( [START_REF] Nguyen | Self-calibration of omnidirectional multi-cameras including synchronization and rolling shutter[END_REF]), and synchronized if needed, using global shutter approximation. This involves a generic Structurefrom-Motion ( [START_REF] Mouragnon | Generic and real-time structure from motion[END_REF]) based on local bundle adjustment, which provides keyframes and matched Harris points that are used latter. We also reduce the accumulated drift of the trajectory thanks to detection and closure of loops inspired by [START_REF] Yousif | An overview to visual odometry and visual slam: applications to mobile robotics[END_REF][START_REF] Younes | Keyframe-based monocular slam: design, survey and future directions[END_REF]: detection using a vocabulary tree of the matched points, closure using global bundle adjustment initialized by pose graph optimizations.

Then the sparse cloud of points is completed for surface reconstruction. For every keyframe, we detect Canny curves and Harris points, match them to those of the previous keyframe ( [START_REF] Lhuillier | Match propagation for image-based modeling rendering[END_REF] and the epipolar constraint. A Canny point is ignored if the angle between its curve tangent and the epipolar line at this point is less than a threshold (π/4). Both Canny and Harris points are reconstructed by ray intersection (using RANSAC and Levenberg-Marquardt methods). We only retain one reconstructed Canny pixel over four consecutive ones in a curve.

Last a pre-filtering is applied to remove gross outliers in the point cloud. Point p i is reconstructed in keyframes located at c j where j ∈ V i (Sec. 2.1). The point is rejected if V i has no more than two elements or all apical angles are small: c j 1 p i c j 2 < where > 0 is a threshold and j k ∈ V i . This roughly means that we reject a 3D point if the ratio of its uncertainty by its depth is greater than a threshold ( [START_REF] Lhuillier | A generic error model and its application to automatic 3d modeling of scenes using a catadioptric camera[END_REF]). The choice of is a trade-off: small values for selecting the scene background and large values for reducing noise in the foreground. We also reject points that are too far in the sky or below the ground as follows. Let v be an estimate of the vertical direction. (It is computed assuming that the multi-camera trajectory is roughly horizontal.) We reject p i if there is j ∈ V i such that altitude v p i is one of the ten smallest or ten largest values in {v p k , j ∈ V k }.

Dataset, parameter setting, and initialization

We take four 1280 × 960 videos at 100Hz by walking during 16 minutes in an university campus using a DIY multi-camera. The multi-camera has four GoPro Hero3 cameras enclosed in a cardboard and mounted on a helmet (top of Fig. 14). The trajectory is about 1500m. long, 3459 keyframes are selected, reconstruction and pre-filtering provide 2.8M points and 13.6M rays. The scene includes streets, cars, facades and vegetation.

We use = π/18 (i.e. 10 • for point selection in Sec. 7.1), α = π/16 (selection of critical edges in Sec. 5.1.2), w 0 = π/2 (P.R. parameter in Sec. 5.1.3). The bottom of Fig. 14 shows the points and the keyframe locations estimated by structurefrom-motion (Sec. 7.1) and the surface using the method of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] (Sec.

Comparison of manifold tests

Here we compare the computation times of the surface reconstruction of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] (summary in Sec. 5.1) by using different manifold tests: the new test based on directed-edges (introduced in Sec. 3) and the old test based on connected components (introduced by [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF]). Both tests check that a vertex is regular in ∂O; their algorithms are detailed in Appendix A.

They obviously provide the same surface. There are three successive operations after Sh.: T.E., C.H.R. and P.R. Using the new test, their times are 3.8s, 27s and 15.6s. Using the old test, their times are 3.9s, 34s and 16.5s. The main acceleration is for the most time expensive operation C.H.R. These tests are done 30M times. Sh. does not change since it has its own test. We also check that the number of edges in D is small (this is useful for the complexity of the new test according to Sec. 3.3): the mean is 8.2, the standard deviation and maximum are 5.9 and 189. The new test is used in the next sections.

Comparison of repairs for singularity removals

Now we compare the method of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] by trying different repairs used by C.H.R.: the new Repair2 based on singularity analysis (introduced in Sec. 4) and the old Repair that adds one tetrahedron at once (introduced by [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF]). We remind that C.H.R. tries to add a lot of sets H ⊆ F \ O to O by using a repair to remove the resulting singular vertices.

Before the use of repair, O/F = 83.54%. Repair2 provides O/F = 85.31% in 10s. Repair provides O/F = 85.21% in 27s. Repair2 has the best ratio and is quite faster than Repair. We note that O/F does not increases a lot, but this is not a reason to omit C.H.R. since it removes visually critical handle. (Fig. 15 shows examples.) Repair2 is used in the next sections.

Comparisons of topologies

Here we experiment three methods using Repair2 and the new manifold test: M0 ( [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF], see also Sec. 5.1), M1 (Sec. 5.2) and M2 (Sec. 5.3). We remind that M1 is M0 with a C.H.R. that can be canceled, M2 is M0 such that T.E. is replaced by C.E.R. and C.H.R. is replaced by S.G. and more selective C.H.R.

Quantitative comparisons

First we check assertions in Sec. 5 about g using Tab. 1. For M0, g is multiplied by about 4.6 due to C.H.R. For M1, g decreases by 35% thanks to C.H.R. For M2, C.E.R. provides g that is about 20% greater than that of T.E. The genus g is the same by S.G. and almost the same by C.H.R. Furthermore, P.R. has similar effects in all cases: g decreases by 10%-16%. At the end, both M1 and M2 provide quite smaller g than M0. (g is divided by 7.6 and 3.7.)

Second we see that there is a trade-off between O/F (which grows if the maximized visibility score function f (O) grows) and small g (i.e. simplified topology): the larger the ratio O/F, the larger the genus g. Indeed, M0 > M2 > M1 for both O/F and g. The trade-off still holds if we use an alternative definition of T.E. by [START_REF] Litvinov | Incremental reconstruction of a complex scene using omnidirectional camera[END_REF], but there are differences: this T.E. provides a quite larger g than the original T.E., which implies that M0 > M1 > M2 for both O/F and g (more details in Appendix F).

Last we check the assumption of Theorem 5 used by M1. It is not 100% meet: the number of connected components of O (which is not c) is 41 after C.H.R. of M1. However, it is 100% meet if we slightly modify C.H.R. (Add condition "v i ∈ c(∂O)" in line 1 of Algorithm 2 of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF].); then we obtain a very similar M1's surface (only 281 different triangles) with same c and g.

Removing holes due to bad or lacking points

First we remind what is a "hole" in our context. Bad input points are remaining after the pre-filtering in Sec. 7.1. Their rays can intersect tetrahedra that become free-space and outside although they should be matter and inside. These tetrahedra form holes. More precisely, there are two kinds of holes: concavities and handles. We explain them for a simple example: a wall in a city such that both planar sides of the wall are seen by the camera trajectory (the c j ). In the first case, a concavity deforms one side of the wall without topology change, e.g. if the wall includes a bad point. Such holes can be removed by P.R. if they are peaks (more details in Appendix C of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF]).

In the second case, a handle connects both wall sides, e.g. if a ray of a bad point crosses the wall.

Second we show examples of holes surrounded by handles (second case) that are removed by our lowered genus constraint. In both Figs. 16 and17, holes surrounded by handles of M0 (column 2) are removed by M1 (column 3) and/or M2 (column 4). M2 usually has the best visual result since M1 sometimes fills large holes that should not be (where there is a lot of free-space F), i.e. its surface connects important scene components that should not be. For example in row 2 of Fig. 17, the space between the tree and the notice board and the ground is in F but the M1 surface immediately connects these three scene components. In row 2 of Fig. 16, we also see that large blunders in F are greatly reduced by the manifold constraint (See M0 result.) and further reduced by the low genus constraint. (See M1 and M2 results.)

Details on C.E.R, S.G. and C.H.R. used in method M2

We note that O/F (Tab. 1) does not increase a lot by operations C.E.R., S.G. and C.H.R. However, this is not a reason to omit them in M2. Fig. 18 shows examples of visual artifacts that occur in the final surface if we omit one of them. If C.E.R. is omitted, a porch forms a blind alley and there is a visual artifact adjacent to a pillar of the porch (row 1). If S.G. is omitted, two adjacent cars are connected by the surface (row 3). If C.H.R. is omitted, spurious handles occur between another post and the ground and a building (row 5). Fig. 18 also shows a top view of all improvements (tetrahedra moved from F \ O to O) done by each of these three operations, which are non negligible.

Unlock shelling blocking

Now we study a method M3, which is M2 with an additional step Unlock (Sec. 6.3) done after C.H.R. and before P.R. Fig. 19 shows examples of visual artifacts of M2 that are removed by M3. It also shows a top view of all improvements (tetrahedra moved from F \O to O) done by Unlock, which are non negligible. The ratio O/F increases from 85.02% to 85.39% thanks to Unlock in 15.3s; the g of M3 is 517. (25% greater than the g of M2.) We also check two assertions in Sec. 6.3. If Unlock does not have the supplementary constraint, there is a lot of topological noise: g of M3 is multiplied by 14 (g = 7279). If we replace the supplementary constraint (before every use of Repair2) by testing δcδg ≥ 0 as in Sec. 5.2.2 (after every successful use of Repair2), g of M3 decreases by 21% (g = 408), the computation time of Unlock is multiplied by 2.1, and O/F = 85.38%.

The video Video.mp4 in the supplementary material shows walkthroughs in the surface generated by M3.

Comparison of accuracies

We compare the accuracies of M0-3's surfaces with respect to a ground truth surface: a synthetic urban scene with textures of real images taken in a city. First a multi-camera video (4 GoPro cameras) is generated by ray-tracing of the scene. The camera trajectory is a 621m long closed loop around buildings. All methods have the same input generated as in Sec. 7.1: 601 keyframes and 687k points. Fig. 20 shows the images at a single location, estimated and ground truth surfaces. A geometric error e between the estimated and ground truth surfaces is computed as follows. Let c i and c g i be the estimated and ground truth locations of the multi-camera at the i-th keyframe. First we set the estimated surface in the coordinate system of the ground truth thanks to the similarity transformation Z that minimizes

E(Z) = 600 i=0 ||Z(c i ) -c g i || 2 . ( 19
)
We obtain √ E(Z)/601 = 9.1cm. Then e(p) is defined by the distance between the point p and the ground truth surface; p is randomly and uniformly sampled in the estimated surface.

Tab. 2 shows c, g and quartiles of e(p). The quartiles below 70% are the sames for all methods. We see that M0-3 have similar accuracies (M2 and M3 are slightly better than M1, which in turn is slightly better than M0.) using a similar number of triangles (1.09M triangles) and times between 36s and 40s. The main difference between methods is the genus g. Its ground truth is 3 since the streets form 3 independent loops and the buildings do not have handles. We see that the topological noise of M0 is quite reduced by M1, M2 and M3. M1 provides the best g but its surface approximates the two streets in the scene center by (completely) blind alleys although they are not. These streets are incompletely blind alleys but have spurious handles in M0, M2 and M3 cases.

Varying α

The angle α is an important parameter since the set of the critical edges depends on α (L α in Eq. 16) and operations C.E.R, C.H.R., S.G. take L α as input. (Sh., T.E. and P.R. do not.) Second Fig. 21 shows surfaces differences for M3 between α = 0.3 and α = 0.1. We see that a large value of α simplifies the topology, e.g. holes/handles are removed in a fence (top). However, this also degrades the O growing at some locations, e.g. near a light post (middle) and below trees (bottom). As in Sec. 7.5.1, there is a trade-off between visibility consistency and topology simplicity.

Comparison with a graph-cut method

First we discuss a graph-cut method (named GC) of [START_REF] Vu | High accuracy and visibility-consistent dense multiview stereo[END_REF]. GC takes ours 3D Delaunay triangulation T as input, then computes a closed surface S ⊂ c(T ) that minimizes a cost function E vis + λE qual composed of a visibility term E vis and 2 a surface quality term E qual , last GC has a Laplacian smoothing. On the one hand, we observe that thin structures (e.g. posts) and foliages are removed from S if the weighting λ increases. On the other hand, S is less robust to "bad" points if λ decreases. These "bad" points can also includes a few good ones which are not enough numerous to generate a descent shape, e.g. a few good points inside a building or a car if the multi-camera is outside. We choose λ = 2 as a trade-off. We also note that a good surface quality in the sense of [START_REF] Vu | High accuracy and visibility-consistent dense multiview stereo[END_REF] implies a low number of singular vertices. Indeed, the percentage of singular vertices of S decreases if λ increases: 1.5% if λ = 1, 0.26% if λ = 2, 0.07% if λ = 3. The computation time of the graph-cut optimization also increases (from 25s to 64s) and the number of triangles in S decreases (from 5M to 3.5M).

Second we compare GC and our method M3 using the default value α = π/16. The number of triangles is 4.2M for M3 and 4.6M for GC. Rows 1-2 in Fig. 22 show that the GC surface is noisier than that of M3 in textured walls. This is coherent with an accuracy evaluation using the dataset in Sec. 7.8: the mean error of M3 is lower than that of GC. (They are 0.33m and 0.43m respectively.) Furthermore, the most complete foliages and trunks of trees are usually obtained by M3 (examples in rows 3-4). However we observe in rows 5-6 that GC is usually the most robust with respect to "bad" points inside buildings. (In row 6, M3 generates a spurious handle connecting two windows at a building corner.) We also note that GC has the best result in row 7 by forcing to O several thin matter tetrahedra that connect roofs of two cars.

Conclusion

Topology constraints (manifoldness, low genus, connectedness) are under-explored in Computer Vision methods that estimate a surface given points reconstructed from images, although these constraints are useful for both surface estimation and applications. This article presents the first surface reconstruction methods that simultaneously enforce visibility consistency and low genus. Starting from a previous work enforcing manifoldness, we quantitatively reduce the genus and show surface improvements including hole removals. A simple modification removes topological noise due to an operation of the original method. A second method is more involved but adjusts the topology simplification thanks to an user parameter during the visibility consistency optimization. Other contributions include an acceleration of the manifold test by using the orientability of the 3D Delaunay triangulation (of the input points), a more efficient removal of surface singularities which improves escapes from local extrema. This removal is based on a study of non-manifoldness near a vertex or an edge of the surface. The local extrema are partly due to limitations of an operation called "shelling" in Combinatorial Topology.

We experiment in a context that we believe useful for applications. First the input point cloud is sparse. This is useful in several contexts including large scale scenes, limited computational resources and initialization of dense stereo. Second the points are reconstructed from videos taken by several consumer cameras (or a spherical camera) mounted on a helmet and by walking (or biking) in complex environments.

Future work includes improvements to escape from local extrema of our optimization problem, e.g. by finding a (provably) set of local operations that reach all manifold surfaces embedded in the 3D Delaunay triangulation. Furthermore, we should improve contour reconstruction, use of connectedness, surface smoothing and visibility score function. Our contributions can also be applied to the incremental surface reconstruction. (We only experiment batch surface reconstruction.) Last a dense stereo refinement should add surface details and reduce reconstruction errors.

Figure 1 :

 1 Figure 1: Examples of manifold surfaces in R 3 having different genuses g. Left: sphere (with dotted equator) meets g = 0. Middle: simple torus meets g = 1. Right: double torus meets g = 2.

Figure 2 :

 2 Figure 2: Local operations in the 2D case. (Tetrahedra are replaced by triangles.) From top to bottom: Shelling (Sh.), Topology Extension (T.E.), Force-Repair (F.R.) and Shrink-Grow (S.G.). Every line shows successive values of O in gray. Sh. adds one tetrahedron at once to O and g does not change. T.E. adds all tetrahedra around a vertex to O so that g can increase. F.R. adds tetrahedra, a singularity is first created and then removed. S.G. first removes then adds tetrahedra to O.

Figure 3 :

 3 Figure 3: Left: triangle orientations induced by the tetrahedron orientation (v 0 , v 1 , v 2 , v 3 ). Front triangles have continuous arrows, back triangles have dotted arrows. Right: two adjacent tetrahedra with consistent orientations (a, b, c, d) and (b, a, c, e). (The triangle abc is duplicated for readability.) The arrows in abc rotate in opposite directions.

  vabc has the following contribution to D: it adds directed edges (b, c), (c, a) and (a, b) to D. If there is no other tetrahedron in O with the vertex v and a triangle face(s) in ∂O, we see that D is a directed cycle and v is regular. Otherwise, D is not a directed cycle (It has supplementary directed edges.) and v is not regular. (There are other triangles t such that v ∈ t ∈ ∂O.) In the second example, tetrahedra vuab, vubc, vuca are in O and triangles vab, vbc, vca are in ∂O. We choose orientations o(vuab) = (v, u, a, b), o(vubc) = (v, u, b, c), o(vuca) = (v, u, c, a) and check orientation consistency, e.g. we have o(vuab) = (u, v, b, a) and o(vubc) = (v, u, b, c) for vuab and vubc. Thus the contribution of the three tetrahedra to D is the following: they add directed edges (a, b), (b, c) and (c, a) to D. Now we do the same observations as in the first example.

Figure 4 :

 4 Figure 4: v ∞ and T ∞ in the 2D case. (Tetrahedra are drawn as triangles, triangles are drawn as edges.) The tetrahedra in T have continuous edges, those in T ∞ \ T have dotted edges. The bold edges form ∂T .

Figure 5 :

 5 Figure 5: T ∞ v \ O (top) and the corresponding g O v restricted to T ∞ v \ O (bottom) in three cases. Case 1 (left): v is regular in ∂O and g O v has one connected component in T ∞ v \ O. Cases 2 and 3 (middle and right): v is singular in ∂O and g O v has two connected component in T ∞ v \ O. The bold edge (middle of top) is singular in ∂O.

Figure 6 :

 6 Figure 6: T ∞ ab ∩ O (top) and the corresponding g O ab (bottom) in three cases. Case 1 (left): ab c(∂O) and g O ab = g ab has one connected component. Case 2 (middle): ab ∈ c(∂O) is not singular in ∂O and g O ab has two connected components. Case 3 (right): ab ∈ c(∂O) is singular in ∂O and g O ab has four connected components.

  abc 2 c 3 ∈ F \ O and abc 3 c 0 ∈ O. Here {abc 2 c 3 } is a connected component of g O ab that can be added to O. If g O a has only one connected component in T ∞ a \ O, this component cannot be added to O since it includes abc 0 c 1 F (and similarly for g O b ). 4.3.2. Algorithm First we need an auxiliary function ReduceSingularity, which is detailed in Algorithm 1. Let V be the set of the singular vertices of ∂O. ReduceSingularity adds to O a tetrahedron set C ⊆ F \ O if this does not add a new vertex to V. In this case, V can decrease and C is added to A. (A is used at the repair end.) Algorithm 1. ReduceSingularity(C,V,O,A) 01: Let {v 0 , v 1 , • • • v c } be the vertex set of C; 02: Let b 0 i be true iff v i V; 03: O ← O ∪ C; 04: Let b 1 i be true iff v i is regular in ∂O; // use manifold test 05: for each i ∈ {0, 1, • • • c} // detect failures 06: if (b 0

Figure 7 :

 7 Figure 7: Repair2 in the 2D case. O is white and F \ O is the set of triangles in gray. From left to right: ∂O is manifold, add H = {abc} to O and obtain the singular vertex set V = {a, b, c}, add to O a connected components of g O a included in F \ O and obtain V = {b}, add to O a connected components of g O b

&& 19 :

 19 ReduceSingularity(C i , V, O, A)) IsImproved= 1; 20: } 21: } while (IsImproved) 22: if (V ∅) { O ← O \ A; return 0; } // restore O if failure 23: // complete O using Shelling if success 24: for each tetrahedron ∆ ∈ A 25: for each 4-neighbor tetrahedron ∆ of ∆ 26: if (∆ ∈ F \ O), try a Shelling started from ∆ ; 27: return 1;

Figure 8 :

 8 Figure 8: Increments (δg, δc) in four cases where we apply O ← O ∪ H. O is white, H is dark gray, T \ O is gray (light and dark). A gray region without (resp. with) a white hole represents a surface in R 3 that has the sphere (resp. torus) topology. The case on the bottom right corner is rejected, the others are accepted.

Figure 9 :

 9 Figure 9: Examples for Theorem 5's proof if k = 2: cases δc = 0 (left) and δc = -1 (right). O is light gray, H is dark gray. In the middle, R 3 \ (|O| \ |∂O|) is white and has three connected components K 0 (unbounded), K 1 and K 2 (bounded). ∂O has three connected components C 0 , C 1 and C 2 . In the left and right, H ⊆ K 2 and ∂(O ∪ H) has at least 2 connected components C 1 and C 0 .

Figure 10 :

 10 Figure 10: Increase (top) or decrease (bottom) of g done by C.E.R. O is gray and F \ O is the set of triangles filled in white. Left: a critical edge e (bold edge) is detected in c(∂O). Middle: a singular vertex (black dot) appears when applying O ← O ∪ (T e ∩ (F \ O)). Right: the singular vertex is removed by repairing (growing) O.

Figure 11 :

 11 Figure 11: Our critical handle removal in the 2D case. O is white and F\O is the set of triangles filled in gray. Left: a critical edge e (bold edges) is detected in c(∂O), π is the dashed line, H is the set of three gray triangles cut by π. Middle: a singular vertex (black dot) appears when applying O ← O ∪ H. Right: the singular vertex is removed by repairing (growing) O.

Figure 12 :

 12 Figure 12: Notations for a shelling blocking. From left to right: O, O , O \ O is homeomorphic to {x ∈ R 3 , ||x|| ≤ 1}, c(O) ∩ c(O \ O) is a disc.

Figure 13 :

 13 Figure 13: Without (left) and with (right) the supplementary constraint in Unlock. We have H = {abcd}, abcd ∈ F \O, F \O forms a coin, O is all around the coin, ab is in the top side and cd is in the bottom side. Left: there are no edges in c(∂O) connecting ab and cd, then Repair2 creates a torus and g increases. Right: bc is in the edge of the coin, then Repair2 does not change g.

  5.1). The surface has 4.2M triangles before sky removal. All experimented methods have the same T , F and initial O by Sh. There are 17.6M tetrahedra in T ; 53% of them are in F. 30% of the vertices in c(∂F) are singular in ∂F. Let O/F be the ratio of the number of tetrahedra in F that are in O. Let c and g be the number of connected components and genus of ∂O. Sh. provides O/F = 83.29%, c = 1 and g = 0 in 13s. (We use a I7-5500U 1600MHz DDR3L laptop.)

Figure 14 :

 14 Figure 14: Top: four GoPro Hero 3 cameras enclosed in a cardboard box and resulting images of our multi-camera in the campus. Middle: top view of the structure-from-motion result including loop closure. Bottom: surface by [22] (best viewed in colors).

Figure 15 :

 15 Figure 15: Examples of removals of spurious handles by [22] without (left) and with (middle and right) C.H.R. involving Repair2.

Figure 16 :

 16 Figure 16: Handles of M0 and their removals by M1, M2 and M3. From left to right: ∂F, M0, M1, M2, M3, and texture of M3. There are views from outside and from inside of a building corner (rows 1&2) and a car (rows 3&4). The arrows are pointing to handles removed by M1/M2/M3.

Figure 17 :

 17 Figure 17: Handles of M0 and their removals by M1, M2 and M3. From left to right: ∂F, M0, M1, M2, M3, and texture of M3. There are views of a small tree (row 1), a notice board under a tree (row 2), an anti-parking post (row 3) and a pillar (row 4). The arrows are pointing to handles removed by M1/M2/M3.

Figure 18 :

 18 Figure 18: Improvements done by C.E.R (top), S.G. (middle) and C.H.R. (bottom) involved in method M2. In each case, there is a top view of tetrahedra switched from F \ O to O (in black) or from O to F \ O (in grey, S.G. only), superimposed by the camera trajectory. There are also local views of the surface obtained without and with the operation.

Figure 19 :

 19 Figure 19: Improvements done by Unlock involved in method M3. Top and middle: local views of the surface without (left) and with (middle and right) Unlock showing three wood posts, two cars, facade and ground, fire hydrant, three and ground. Bottom: a top view of tetrahedra switched from F \ O to O (in black) by Unlock, superimposed by the camera trajectory.

Figure 20 :

 20 Figure 20: Experiments using a synthetic urban scene. Top: images at a location. Middle: top views of SfM result and M3 surface. Bottom: local views of ground truth and M3 surfaces.

Figure 21 :

 21 Figure 21: Varying α for the method M3: α = 0.3 on the left and α = 0.1 on the middle and right. There are views of a fence in front of a building, a light post between two streets, and tall grass under trees.

First

  we detail a link between g and α. By varying α around its default value π/16 ≈ 0.2, Tab. 3 shows that g increases if α decreases for M0, M2 and M3. Indeed, new handles can be generated at the neighborhood of every critical edge in L α where operations C.E.R./C.H.R. are applied and the number of edges in L α increases if α decreases. In contrast to this, g slightly decreases if α decreases for M1. There is a reason: the cancellations of F.R. involved in C.H.R. of M1 prohibit g increments. (This is the only differences between M0 and M1.) In all cases, if α decreases from 0.3 to 0.1., c increases in range[START_REF] Ziegler | Shelling polyhedral 3-balls and 4-polytopes[END_REF] 70] and the computation time of C.H.R. is multiplied by about 3.5.

Figure 22 :

 22 Figure 22: Differences between GC (left) and M3 (middle and right). M3 is the best in rows 1-4 and GC is the best in rows 5-7 (details in the text). All views are from the outside except those in rows 2 and 5.

Figure B. 23 :

 23 Figure B.23: Notations for Theorem 2 proof with l = 1 and m = 4. Left: the triangles in ∂O including the vertex v. Here v is regular in ∂O. Right: the tetrahedra in T uv ∩ O. We have uva i a i+1 ∈ O where 1 ≤ i < 4, vua 1 ∈ ∂O, vua 4 ∈ ∂O. In both cases, the bold edges are the U edges, i.e. the v-opposite edges in triangles of ∂O.

  and others. (This is not the topic of the paper.) Intuitively, a large O in F in the inclusion sense is a good solution if ∂O is manifold. Every operation tries to add to O tetrahedra that are in F \ O. It iteratively generates a series of tetrahedron sets O 0 , O 1 • • • O n where O 0 is the initial value of O and n > 0. The operation is successful iff ∂O n is manifold and f (O n ) ≥ f (O 0 ). It returns the new value of O: O n if success or O 0 if failure.

	2.2.2. Operations Sh., T.E., F.R., C.H.R.s, and S.G.
	Almost all operations are local. Fig. 2 summarizes them in
	the 2D case: tetrahedra are replaced by triangles. (The latter can

  5.2.1. Efficient computation of δvδe + δt Let H ⊆ F \ O be all tetrahedra that are added to the initial O by a successful F.R. operation during C.H.R. (Note that H ⊆ H by using the notation H in Sec. 2.2.2.) Let

Table 1 :

 1 Labeling, surface topology (number of connected components and genus) and times of methods M0-2. Reminders: both M0 and M1 successively apply Sh.-T.E.-C.H.R.-P.R., M2 successively applies Sh.-C.E.R.-S.G.-C.H.R.-P.R.

	step	method	O/F	c	g	time
	T.E.	M0	83.54% 368 377	3.9s
	C.H.R.	M0	85.31% 483 1759 10.4s
	P.R.	M0	79.87% 56 1549 15.6s
	T.E.	M1	83.54% 368 377	3.9s
	C.H.R.	M1	84.75% 499 243 16.5s
	P.R.	M1	79.14% 56	204	16s
	C.E.R.	M2	84.90% 154 454	4.9s
	S.G.	M2	84.99% 155 454	5.1s
	C.H.R.	M2	85.02% 171 459	5.6s
	P.R.	M2	79.45% 55	411 15.6s

Table 2 :

 2 Mean, standard deviation and quartiles of geometric error e(p) in cm, c and g for the methods M0-3.

	M.	ē	σ e	70% 80% 90% c g
	M0 34.73	103	14	26	76	3 82
	M1 34.12	101	14	26	78	3 2
	M2 32.89 100.8	14	24	70	3 17
	M3 32.93 100.8	14	24	70	3 33

Table 3 :

 3 Genuses obtained by the methods M0-3 and values of α.

	α	0.1	0.15	0.2	0.25 0.3
	M0 3774 2269 1531 1156 945
	M1 194	199	207	220 225
	M2 920	580	406	308 237
	M3 1019 669	511	417 353

using equalities (a, b, c, d) = (d, b, a, c) and (b, a, c, e) = (e, a, b, c).

Addendum: using α vis = 1.

Acknowledgments

Thanks to CNRS, Université Clermont Auvergne and Institut Pascal for funding Maxime Lhuillier. Thanks to the CRISTAL project for providing the ground truth VRML model used in the quantitative evaluation.

Appendix A. Algorithms of manifold tests

Here we use shortened notations:

In our implementation, every tetrahedron ∆ ∈ T ∞ (Sec. 4.1) stores four booleans t ∆ i ∈ ∂O in addition to a boolean ∆ ∈ O. This redundancy accelerates all tests by reducing the use of indices or pointers to neighboring tetrahedra.

Appendix A.1. Our manifold test

We remind that π ∈ A 4 iff (π0, π1, π2, π3) ∈ Ã4 where Ã4 = { (0, 1, 2, 3), (0, 2, 3, 1), (0, 3, 1, 2),

(1, 2, 0, 3), (1, 0, 3, 2), (1, 3, 2, 0), (2, 3, 0, 1), (2, 0, 1, 3), (2, 1, 3, 0), (3, 0, 2, 1), (3, 2, 1, 0), (3, 1, 0, 2) }.

(A.1) Algorithm 3 presents (in C style) our directed edge-based manifold test (Sec. 3.3). In the first step (collect the directed edges), lines 3-4 select ∆ ∈ T v ∩ O, then lines 5-20 describe every π ∈ A 4 as follows. Lines {5,9,13,17} select the four cases of π0 ∈ {0, 1, 2, 3}. Lines 6-8 are the three sub-cases of π0 = 0 based on tests t ∆ π1 ∈ ∂O where π1 ∈ {1, 2, 3}, then (v ∆ π2 , v ∆ π3 ) is collected in every sub-case such that π is even. The other cases of π0 are similar. In the second step (Check that D is a directed cycle.), we permute the edges in E to check that D is a closed path, then check that every vertex is traversed only once by this path. [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF]), that we experiment and compare to Algorithm 3. It checks that the vertex v is regular in ∂O by a traversal of the adjacency graph of T ∞ v (proof in Theorem 3 of [START_REF] Lhuillier | 2-manifold tests for 3d delaunay triangulation-based surface reconstruction[END_REF]). Every tetrahedron ∆ has indices (or pointers) to its four adjacent tetrahedra

Tetrahedron-based manifold test for v 01: // every tetrahedron ∆ has boolean f ∆ = 0 02: for each tetrahedron ∆ ∈ T ∞ v 03: f ∆ = 1; 04: c= 0; 05: let P be a table of tetrahedra indices, and let n= 0; 06:

First we show Lemma 1.

Proof. First assume that edge ab ∈ U and show that triangle abv ∈ ∂O. Thus (a, b) ∈ D or (b, a) ∈ D. In both cases, Eq. 8 implies that abv ∈ ∂O.

Second assume that abv ∈ ∂O and show that ab ∈ U. There is a tetrahedron ∆ ∈ O such that abv is a face of ∆. Let x be a vertex such that ∆ = vxab. Since there are exactly two possible orientations, Proof. According to [START_REF] Boissonnat | Geometric structures for three-dimensional shape representation[END_REF][START_REF] Lhuillier | 2-manifold tests for 3d delaunay triangulation-based surface reconstruction[END_REF], a vertex v ∈ c(∂O) is regular in ∂O iff the v-opposite edges in the triangles of ∂O (having v as vertex) form a cycle. Now we use Lemma 1 and obtain the result.

We also need the two following lemma.

Assume (reductio ad absurdum) that the former is true. Thanks to Theorem 1, ∆ and ∆ are consistently oriented by using orientations o(∆) and o(∆ ): there are π ∈ A 4 and π

Then we show that

). Now we use Eq. B.1 with adjacent tetrahedra uva i+1 a i and uva i+1 a i+2 , which are consistently oriented, and obtain (v, u, a i+1 , a i+2 ) = o(vua i+1 a i+2 ).

Last we show Theorem 2.

Proof. Assume that D is a directed cycle and show that v is regular in ∂O. U is a cycle since D is a directed cycle. Then Lemma 2 implies that v is regular in ∂O.

Conversely, assume that v is regular in ∂O and show that D is a directed cycle. Thanks to Lemma 2, U is a cycle. Since D cannot include two inverse edges (Lemma 6), we only need to show that there are edges (a, u) and (u, b) in D for every vertex u ∈ c(U).

Let a vertex u ∈ c(U). Since T is a 3D Delaunay triangulation, there are distinct vertices a i such that

with a possible exception a 0 = a k+1 . Thanks to Lemma 1, if uva i a i+1 ∈ O and uva j ∈ ∂O and j ∈ {i, i + 1}, ua j ∈ U. Since U is a cycle, the number of such uva j is exactly 2. Thus we have l and m such that

and vua l ∈ ∂O and vua m ∈ ∂O. ) is (v, u, a l , a l+1 ) or (u, v, a l , a l+1 ).

In case 1, (v, u, a l , a l+1 ) = o(vua l a l+1 ) and we show that (u, a l ) ∈ D and (a m , u) ∈ D. According to Eq. 8 and since (v, a l+1 , u, a l , ) = (v, u, a l , a l+1 ) and vua l ∈ ∂O, (u, a l ) ∈ D. Thanks to Lemma 7 and (v, u, a l , a l+1 ) = o(vua l a l+1 ), we have (v, u, a m-1 , a m ) = o(vua m-1 a m ). Since (v, a m-1 , a m , u) = (v, u, a m-1 , a m ) and vua m ∈ ∂O, Eq. 8 implies (a m , u) ∈ D.

In case 2, (u, v, a l , a l+1 ) = o(vua l a l+1 ) and we show that (a l , u) ∈ D and (u, a m ) ∈ D. According to Eq. 8 and since (v, a l+1 , a l , u) = (u, v, a l , a l+1 ) and vua l ∈ ∂O, (a l , u) ∈ D. Thanks to Lemma 7 and (u, v, a l , a l+1 ) = o(vua l a l+1 ), we have

Appendix D. Shrink-Grow (S.G.) algorithm Algorithm 5 is Shrink-Grow. (S.G. is summarized in Sec. 2.2.2.) This algorithm is repeated until f (O) does not change or a maximum number of iterations is reached. It uses a function Growing(∆) which updates O by Shelling using ∆ as the first tetrahedron that we try to add to O (i.e. Algorithm 1 of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] with inputs Q 0 = {∆} and O ∅). The function Growing also returns the set of tetrahedra that it adds to O. Algorithm 5. Shrink-Grow (one iteration) 01: Let L α be defined by Eq. 16; 02:

for each tetrahedron ∆ ∈ L seed 11:

L add ← L add ∪ Growing(∆); 12:

R sub = ∆∈L sub f (∆); // f is defined in Sec. 2.2.1 13:

R add = ∆∈L add f (∆); Let P ⊆ F be defined (as in Sec. 6.1 of [START_REF] Lhuillier | Overview of shelling for 2-manifold surface reconstruction based on 3d delaunay triangulation[END_REF]) by

using distinct vertices a 0 , a 1 , b 0 , b 1 , c 0 and c 1 . We summarize P: it has one "internal" tetrahedron a 0 a 1 b 1 c 0 since every triangle face of a 0 a 1 b 1 c 0 is included in another tetrahedron of P having vertex b 0 or c 1 . Furthermore,

be the set of the singular vertices of ∂O i .

21 Appendix F. Use T.E. of [START_REF] Litvinov | Incremental reconstruction of a complex scene using omnidirectional camera[END_REF] instead of T.E. of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF] In this paper, we use T.E. of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse structure-from-motion data[END_REF]: if v ∈ c(∂O) and T v ⊆ F and ∂(O ∪ T v ) is manifold, then update O ← O ∪ T v . However there is another T.E. in algorithm IV.2 of [START_REF] Litvinov | Incremental reconstruction of a complex scene using omnidirectional camera[END_REF]: if v ∈ c(∂O) and ∂(O ∪ (T v ∩ F)) is manifold, then update O ← O ∪ (T v ∩ F). Since the latter is less constrained than the former, it provides the largest O growing. (See Tab. F.4.) However the latter generates a lot of topological noise (g is multiplied by 5.7.) and its computation time is 5 times larger than that of the former. As a consequence, the g of M0 increases by 72% and the g of M1 is multiplied by 4. M2 is unchanged since it does not use T.E. As in Sec. 7.5.1, we see a trade-off between visibility consistency and topology simplicity: M0 > M1 > M2 for both O/F and g.

According to