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

Tests for 3D Delaunay Triangulation-based Surface Reconstruction

Keywords: Reconstruction, Volumetric models, Geometric topology, Duality and planar graphs

This is a companion paper of a previous work on the surface reconstruction from a sparse cloud of points, which are estimated by Structure-from-Motion. The surface is a 2-manifold sub-complex of the 3D Delaunay triangulation of the points. It is computed as the boundary of a list of tetrahedra, which grows in the set of Delaunay tetrahedra. Here we detail the proofs for the 2-manifold tests that are used during the growing: we show that the tetrahedron-based test and the test for adding (or subtracting) one tetrahedron to (or from) the list are equivalent to standard tests based on triangles.

Introduction

Recently, batch [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF] and incremental [START_REF] Yu | Incremental reconstruction of manifold surface from sparse visual mapping[END_REF][START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional structure-from-motion data[END_REF] methods have been introduced to reconstruct a surface from sparse Structure-from-Motion points. These methods adaptively divide up the space using a 3D Delaunay triangulation (tetrahedralization) of the points, whose tetrahedra are labeled as free-space or matter using visibility constraints provided by Structure-from-Motion. The surface is then computed as the boundary of a list O of tetrahedra, which grows in the set of free-space tetrahedra. This list must be as large as possible subject to the constraint that the boundary is a 2-manifold.

These methods represent the 3D Delaunay triangulation by the adjacency graph g of the tetrahedra, since this is adequate for both memory space [START_REF] Boissonnat | Triangulations in CGAL[END_REF] and tetrahedron labeling. The surface is implicitly represented by the set which includes every g edge linking a tetrahedron in O and another one not in O. However, the standard 2-manifold tests for a surface are defined on a surface-based representation (e.g. the adjacency graph of triangles), not a volumetric one as in ours.

Thus, two 2-manifold tests directly based on tetrahedra are used in [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF] for efficiency. The former checks that the surface is a 2-manifold near surface vertex v by counting the connected components of a g sub-graph of the tetrahedra having v as vertex. The latter checks that a single tetrahedron can be added to O without loss of the 2-manifold property.

However, the complete proofs of the two tests have not been published before. The former has a short and sketch proof in Appendix A of [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF], which ignores the case of a surface with a vertex at the boundary of the Delaunay triangulation. The latter only has a visual proof [START_REF] Boissonnat | Geometric structures for threedimensional shape representation[END_REF] (a figure) or a partial proof [START_REF] Gezahegne | Surface reconstruction with constrained sculpting[END_REF] (sufficient but not necessary condition) in the special case where the surface encloses all input points. Here we give detailed proofs with necessary and sufficient conditions in all cases. The paper provides prerequisites (Sec. 2), statements (Sec. 3) and proofs (Sec. 5 and 4) of the tests.

Prerequisites (Main Definitions, Properties)

The majority of prerequisites are in [START_REF]Handbook of discrete and computational geometry[END_REF][START_REF] Giblin | Graphs, surfaces and homology[END_REF]. We introduce integers k ≥ 0 and n > 0. A (geometric) simplex σ is the convex hull of k + 1 points v 0 • • • v k in general position in R n , i.e. v 1v 0 • • • v kv 0 are linearly independent. We say that σ ′ is a face of σ if σ ′ is the convex hull of some of the v i above (thus σ ′ ⊆ σ). A simplicial complex K in R n is a finite set of simplices in R n such that 1. if σ ∈ K and σ ′ is a face of σ, σ ′ ∈ K. 2. if {σ, σ ′ } ⊆ K, σ ∩ σ ′ is empty or a face of σ and σ ′ . Let V ′ be a set. An abstract simplicial complex S is a finite set of subsets of V ′ such that A ∈ S and B ⊂ A imply B ∈ S. It is implicitly defined by every simplicial complex K as follows: V ′ is the vertex set of K and S is the family of the vertex sets of the K simplices. Conversely, we say that K is a realization of S in R n . In this paper, we often use the same notation for a simplicial complex and its abstract version.

The elements of S are called (abstract) simplices. Their faces are their subsets. If simplex σ has exactly k+1 vertices, σ has a dimension k. Simplex σ is a vertex, edge, triangle or tetrahedron if k is 0,1,2 or 3, respectively. We use bold fonts for vertices in an (abstract) simplicial complex, e.g. a is a vertex, ab is an edge, abc is a triangle. We say that edge ab is c-opposite in triangle abc, triangle

abc is d-opposite in tetrahedron abcd, σ is σ ′ -incident if σ ′ ⊂ σ.
Let L ′ be a subset of simplices in an (abstract) simplicial complex. The closure c(L ′ ) is the set of the faces of the L ′ simplices; c(L ′ ) is an (abstract) simplicial complex which includes L ′ . Let σ ∈ L ′ . If σ is not included in another simplex in L ′ , we say that σ is maximal in L ′ . In a k-dimensional pure (abstract) simplicial complex, every maximal simplex has dimension k, therefore every simplex is included in a simplex of dimension k.

The boundary ∂σ of a k-dimensional simplex σ is the set of its k+1 faces of dimension k-1: a tetrahedron has four triangles, a triangle has three edges, an edge has two vertices, a vertex has an empty boundary. If all simplices in L ′ are tetrahedra, boundary ∂L ′ is the set of triangles such that every triangle is a face of exactly one tetrahedron in L ′ . In this case, c(L ′ ) and c(∂L ′ ) are 3D and 2D pure (abstract) simplicial complexes, respectively.

Let V be a set of m ≥ 4 points in R 3 . Let T be a set of tetrahedra which meets three conditions:

1. V is the vertex set of the tetrahedra in T . 2. the convex hull of V is the union of tetrahedra in T . 3. the circumscribing sphere of every tetrahedron in T does not contain a V vertex in its interior.

Thus c(T ) is a 3D pure simplicial complex in R 3 , which is called a 3D Delaunay triangulation of V . Every triangle in c(T ) is included in exactly two tetrahedra in c(T ), except those in ∂T (we have ∂T = ∅).

To avoid special cases in the statements and proofs of our tests, we extend c(T ) to a 3D pure abstract simplicial complex where every triangle is included in exactly two tetrahedra [START_REF] Boissonnat | Triangulations in CGAL[END_REF].

Its vertex set is V ∞ = V ∪ {v ∞ }, where new vertex v ∞ is called the infinite vertex (v ∞ / ∈ R 3 ). Its tetrahedron set is T ∞ = T ∪ {abcv ∞ , abc ∈ ∂T } (here we use the abstract versions of c(T ) simplices). Thus c(T ∞ ) is a 3D pure abstract simplicial complex including the abstract version of c(T ). Assume σ ∈ c(T ∞ ). We say that σ is infinite if v ∞ is a vertex of σ, otherwise σ is finite. Let k ∈ {2, 3} and L be a subset of k-dimensional simplices in c(T ∞ ). If k = 3, we use notation L c = T ∞ \ L.
Let τ be a vertex or an edge in c(T ∞ ). We also use notation

L τ = {σ ∈ L, τ ⊂ σ}. For example, T ∞ τ is the set of the τ -incident tetrahedra in T ∞ , L c
τ is the set of the τ -incident tetrahedra in T ∞ \L. If τ is the edge with vertices a and b, L τ = L a ∩ L b . Indeed, a triangle (or a tetrahedron) includes edge ab iff it includes vertices a and b ("iff" means "if and only if"). Assuming vertex v ∈ c(T ), v ∈ c(∂T ) iff T ∞ v contains at least one infinite tetrahedron. We also have ∂T ∞ v = {abc, abcv ∈ T ∞ }. A cycle is a graph with m vertices q 1 • • • q m and m edges q 1 q 2 , • • • , q m-1 q m , q m q 1 where m > 2. Let (V, E) be the graph with vertex set V and edge set E. Let {V 1 , V 2 } be a partition of V . Let F be the set of all edges in E such that every edge in F has a vertex in

V 1 and a vertex in V 2 . Set F is a cut in (V, E). Graph (V, E \ F ) has at least two connected components: one in V 1 and another one in V 2 . Set F is a minimal cut [3] if V 1 and V 2 are two connected components in (V, E \ F ).
Two simplices σ 1 and σ 2 in L are adjacent if σ 1 ∩ σ 2 has dimension k -1 (e.g. two tetrahedra are adjacent if their intersection is a triangle). Let g L be the adjacency graph of L: the vertices of g L are the simplices in L, every edge of g L links two distinct vertices in g L if the two corresponding simplices in L are adjacent. Let τ be a vertex or and edge in c(T ), and g be the adjacency graph of T ∞ . The adjacency graph of T ∞ τ is shortened g τ ; g τ is a connected sub-graph of g. If τ is an edge, g τ is a cycle. If vertex v is in c(T ), the vertices and edges of c(∂T ∞ v ) form a connected graph. Proofs for these properties are in the supplementary material (available at http://maxime.lhuillier.free.fr/pJmiv15Supp.pdf).

A n-ball B has center x ∈ R n and radius r > 0, i.e.

B = {y ∈ R n , ||x -y|| < r} where ||.|| is the Euclidean norm. B is open in R n . If L ⊆ c(T ), we define |L| = ∪ σ∈L σ (union of convex hulls). A set included in R n (e.g. |L|)
has the topology induced by the n-balls [START_REF] Zomorodian | Topology for computing[END_REF]. Topological spaces X and Y are homeomorphic if there is a bijective and continuous function ϕ such that ϕ -1 is continuous and ϕ

(X) = Y ; ϕ is a homeomorphism. Let M ⊆ R 3 . We say that x ∈ M is regular in M if x has a neighborhood in M homeomorphic to a 2-ball. Intuitively, M has a local 2D parametrization at x. If every x ∈ M is regular in M , M is a 2-manifold in R 3 .

Overview of 2-Manifold Tests

Let O be a list of tetrahedra included in T ∞ such that every triangle in ∂O is finite. We would then like to check that |∂O| is a 2-manifold. We use two definitions.

Definition 1 (Good Edges) An edge in c(∂O) is ∂O- good if it is included in exactly two triangles of ∂O. Definition 2 (Good Vertex) A vertex v in c(∂O) is ∂O-good if the v-incident triangles in ∂O can be or- dered as t 0 , t 1 , • • • t k such that t i ∩ t (i+1) mod (k+1) is an edge for every i ∈ {0, 1, • • • k}.
According to [START_REF] Vegter | Computational Topology[END_REF],

Theorem 1 (Global Test) |∂O| is a 2-manifold iff the vertices and the edges in c(∂O) are ∂O-good.

We have the same test for a single vertex: For the completeness of the paper, we also provide proofs of these theorems in the supplementary material. Fig. 1 shows two examples where O has two tetrahedra and both conditions (good vertices and good edges) are needed to obtain a 2-manifold. On the left, vertex v is not regular in |∂O| since v is not ∂O-good. On the right, every point of edge e is not regular in |∂O| since e is not ∂O-good. All other vertices and edges are ∂O-good.

Theorem 2 (Triangle-based Test) A vertex v in c(∂O) is regular in |∂O| iff v and the v-incident edges in c(∂O) are ∂O-good.
Since Theorem 1 condition is the conjunction of the conditions of Theorem 2 for all vertices in c(∂O),

Corollary 1 |∂O| is a 2-manifold iff every vertex of c(∂O) is regular in |∂O|.
The Triangle-based Test can be rewritten using edges [START_REF] Boissonnat | Geometric structures for threedimensional shape representation[END_REF]:

Corollary 2 (Edge-based Test) A vertex v in c(∂O)
is regular in |∂O| iff the v-opposite edges in the triangles of ∂O having v as vertex form a cycle. Now assume that our data structure is the adjacency graph g of the tetrahedra of T ∞ . Boundary ∂O is represented by a cut: the g edges between O and O c . In this case, the following test is preferred [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF] over the Triangle-based Test.

Theorem 3 (Tetrahedron-based Test) Let g O v be the graph obtained from g v by removing the edges between a tetrahedron in O and another in

O c . Vertex v ∈ c(∂O) is regular in |∂O| iff g O v exactly has 2 con- nected components. These components are O v and O c
v . The proof is in Sec. 5. In [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional structure-from-motion data[END_REF][START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF], the implementation of this test is a graph traversal of g O v . Fig. 2 shows examples to experiment the above tests: if f = 0 and every

1. v / ∈ c(∂O), g O v is connected since g O v = g v . 2. v is not regular in |∂O|, the v-opposite edges does not form a cycle, g O v has 3 connected components. 3. v is regular in |∂O|, the v-opposite edges form a cycle, g O v has
: v / ∈ c(∂O) and O v = T ∞ v . Middle: v is non regular and O v = T ∞ v \ {abcv, cdev}. Bottom: v is regular and O v = T ∞ v \ {abcv}. Left: tetrahedra in O v . Middle-left: v-incident triangles in ∂O. Middle-right: graph of the v-opposite edges in triangles of ∂O. Right: graph g O v (a circle is filled iff the corresponding tetrahedron is in O c v ). f=4 f=0 f=1 f=2 f=3
∆ vertex v meets O v = ∅. -if f = 1 and the ∆ vertex v, which is opposite to the ∆ triangle included in ∂O, meets O v = ∅. -if f = 2 and the ∆ edge e, which is not in one of the two ∆ triangles included in ∂O, meets O e = ∅. -if f = 3. -if f = 4.
Fig. 3 shows the five cases. The proof is in Sec. [START_REF] Gezahegne | Surface reconstruction with constrained sculpting[END_REF]. In [START_REF] Yu | Incremental reconstruction of manifold surface from sparse visual mapping[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 and omnidirectional structure-from-motion data[END_REF], we check that O v = ∅ or O e = ∅ by a single reading of T ∞ v , the set of v-incident tetrahedra (for edge e with vertices v and w, we look for the tetrahedra in T ∞ v having vertex w).

Since ∂O = ∂O c and subtracting one tetrahedron from O is like adding one tetrahedron to O c , we can rewrite the addition Theorem 4 as a subtraction one. This result is used in [START_REF] Litvinov | Incremental solid modeling from sparse and omnidirectional structure-from-motion data[END_REF] (without proof). A similar result is also used in [START_REF] Boissonnat | Geometric structures for threedimensional shape representation[END_REF][START_REF] Gezahegne | Surface reconstruction with constrained sculpting[END_REF] where only cases f ∈ {1, 2} occur.

4 Proof for "Adding One Tetrahedron" First we study the simplicial complex "between O and ∆" in the two following lemmas. Intuitively, it is indifferently defined by the intersection of closures of O and ∆, or by the intersection of closures of ∂O and ∂∆. We have O ⊂ T ∞ , ∂O ⊂ c(T ) and ∆ ∈ T \ O in the paper.

Lemma 1 c(O) ∩ c(∆) is a simplicial complex in R 3 . Furthermore, c(O) ∩ c(∆) = c(∂O) ∩ c(∂∆). As a con- sequence, the triangles in c(O) ∩ c(∆) are exactly those in ∂∆ ∩ ∂O. Lemma 2 |∂∆|∩|∂O| = |c(∂O)∩c(∂∆)|. If c(O)∩c(∆) is 2D pure, |∂∆| ∩ |∂O| = |∂∆ ∩ ∂O|.
We also examine the triangles in ∂(O ∪ {∆}):

Lemma 3 ∂(O ∪ {∆}) = (∂O ∪ ∂∆) \ (∂O ∩ ∂∆).
The proofs of the three lemmas above are in Appendix A. They can be omitted at the first reading since they are essentially consequences of basic properties of (abstract) simplicial complexes.

In Lemma 4, we convert the condition of Theorem 4 to a more tractable condition for our proof. It is also useful to visualize all cases in Theorem 4 at once.

v of ∆ such that O v = ∅, v ∈ c(O) ∩ c(∆).
There is therefore a maximal simplex τ in c(O) ∩ c(∆) which contains v, and τ is not a triangle. Conversely, let τ be a maximal

simplex in c(O) ∩ c(∆) (which is not a triangle). Every vertex v of τ is in simplicial complex c(O) ∩ c(∆). This implies that v is a vertex of ∆ and O v = ∅. Assume f = 1. We have ∆ = abcv and ∂O ∩ ∂∆ = {abc}. According to Lemma 1, abc is the unique tri- angle in c(O) ∩ c(∆). If O v = ∅, v ∈ c(O) ∩ c(∆).
There is a maximal simplex τ in c(O) ∩ c(∆) which contains v, and τ is not a triangle (otherwise τ = abc, which is impossible). Conversely, let τ be a maximal simplex in c(O)∩c(∆) which is not a triangle. Assume (reductio ad absurdum) that v is not a vertex of τ . Since τ ⊂ abcv, τ ⊆ abc. Since τ is maximal, τ = abc (impossible).

|∂(O ∪{∆})| is not a 2-manifold since τ is not ∂(O ∪ {∆})-good. First we detail (∂(O ∪ {∆})) τ thanks to Lemma 3 and (∂O ∪ ∂∆) \ (∂O ∩ ∂∆) = (∂O \ ∂∆) ∪ (∂∆ \ ∂O). Using shortened notations L O τ = (∂O \ ∂∆) τ and L ∆ τ = (∂∆ \ ∂O) τ , we have (∂(O ∪ {∆})) τ = L O τ ∪ L ∆ τ . Assume that τ is a vertex. Since τ ∈ c(∂O) ∩ c(∂∆) (Lemma 1), there are triangles t O and t ∆ such that τ ⊂ t O ∈ ∂O and τ ⊂ t ∆ ∈ ∂∆. Furthermore t O / ∈ ∂∆ (otherwise, τ ⊂ t O ∈ ∂O ∩ ∂∆ ⊆ c(O) ∩ c(∆)
and thus τ is not maximal) and similarly t ∆ / ∈ ∂O. We obtain

L ∆ τ = ∅ and L O τ = ∅. Let t ′ ∆ ∈ L ∆ τ , t ′ O ∈ L O τ and e = t ′ O ∩ t ′ ∆ .
If e is an edge, τ ⊂ e ∈ c(O) ∩ c(∆) and τ is not maximal (impossible). Thus the adjacency graph of

L ∆ τ ∪ L O τ is not connected: τ is not ∂(O ∪ {∆})-good. If τ is an edge, g τ is a cycle of adjacent tetrahedra ∆ 0 ∆ 1 • • • ∆ n ∆ 0 in T ∞ τ (n ≥ 2) where ∆ 0 = ∆. We have ∆ 1 / ∈ O and ∆ n / ∈ O (otherwise there is i ∈ {1, n} such that τ ⊂ ∆ i ∩ ∆ ∈ c(O) ∩ c(∆) and τ is not maximal).
Since O τ is included in this cycle and ∆ / ∈ O τ = ∅, there are (at least) four pairs (∆ i , ∆ (i+1) mod (n+1) ) such that one tetrahedron is in O ∪ {∆} and the other one is not. Thus, four triangles in

∂(O ∪ {∆}) include τ : τ is not ∂(O ∪ {∆})-good. ⊓ ⊔ Lemma 6 If |∂O| is a 2-manifold and c(O) ∩ c(∆) is a 2D pure simplicial complex, |∂(O ∪ {∆})| is a 2- manifold.
Proof Let v be a vertex in c(∂O). Let G be the graph with the vertices V and the edges E in c(∂T ∞ v ). Let V * be the triangles in c(∂T ∞ v ), i.e. V * = ∂T ∞ v . Let G * be the adjacency graph of the triangles in V * . Let E * be the G * edges, i.e. e * ∈ E * is an edge between triangles {t j , t k } ⊆ V * if t j ∩ t k is an edge e ∈ E. The top of Fig. 4 shows G and G * if T ∞ v is a set of 6 tetrahedra. Here we need an additional definition. A drawing of graph G = (V, E) is a function ϕ defined on V ∪ E such that V and ϕ(V ) ⊆ R 2 are bijective, ϕ(ab) is a non self-intersecting curve in R 2 whose endpoints are ϕ(a) and ϕ(b) if ab ∈ E. Furthermore, two such curves cannot intersect except at their endpoints.

The proof of Theorem 3 has six steps:

1. show that e → e * is a bijection between E and E * .

2. find a G drawing using a realization of c(T ∞ v ) in R 3 . 3. show that G * has a drawing which is dual [START_REF] Diestel | Graph theory[END_REF] to G drawing (Fig. 4 shows one example for G and G * ). 4. express Theorem 3 in terms of G * and a minimal cut of G * . 5. express Corollary 2 in terms of G and a cycle of G. 6. show that Theorem 3 and Corollary 2 are equivalent through duality [START_REF] Diestel | Graph theory[END_REF] between cycles of G and minimal cuts of G * .

Step 1. We have

∂T ∞ v = {abc, abcv ∈ T ∞ v }.
Therefore g v (the adjacency graph of tetrahedra in T ∞ v ) and G * (the adjacency graph of triangles in ∂T ∞ v ) are the same. Since every triangle abv is included in exactly two tetrahedra in T ∞ , every edge ab ∈ E is included in exactly two triangles in V * . The function e → e * is thus well defined and bijective between E and E * .

Step 2. First we need a realization

K of c(T ∞ v ) in R 3 . If v / ∈ c(∂T ), T ∞ v does not contain infinite tetra- hedra and c(T ∞ v ) is a simplicial complex. We use K = c(T ∞ v ). If v ∈ c(∂T ), T ∞ v contains infinite tetrahedra. In this case, another K is required. Lemma 7 c(T ∞ v ) has realization K in R 3 if v ∈ c(∂T ).
Proof Note that c(T ∞ v ) is the union of two abstract simplicial complexes: c(T v ) and the closure of the vv ∞incident tetrahedra in c(T ∞ ). The former is also a simplicial complex in R 3 . The idea of the proof is the following. We find w ∈ R 3 such that the latter has realization K ′ in R 3 by replacing v ∞ by w, and such that K ′ ∪c(T v ) is a simplicial complex. We obtain K = K ′ ∪ c(T v ).

Let J = {vab, vab ∈ ∂T } and w ∈ R 3 . Assume 1. if vab ∈ J, {w, v, a, b} are in general position.

Let w * J = {wσ, σ ∈ c(J)} where wσ is the geometric simplex whose vertex set is that of σ plus w. Assume 2. if {wσ, wσ ′ } ⊆ w * J, (wσ) ∩ (wσ ′ ) is empty or a face of wσ and wσ ′ .

Thus K ′ = {w} ∪ c(J) ∪ (w * J) is a simplicial complex: it is the cone [START_REF] Giblin | Graphs, surfaces and homology[END_REF] on c(J) with vertex w. Assume We have (1) since w / ∈ ∪ t∈J π t ; ( 2) is a consequence of (wσ)∩(wσ ′ ) = w(σ∩σ ′ ). We have w(σ∩σ ′ ) ⊆ (wσ)∩ (wσ ′ ). Conversely, let x be a point in (wσ) ∩ (wσ ′ ) and show x ∈ w(σ ∩ σ ′ ). Let a be a point in σ and b be a point in σ ′ such that x ∈ wa∩wb. Assume (reductio ad absurdum) that a = b. Since w, x, a, b are collinear, a ∈ wb. Let t be a triangle in J such that b ∈ t. Thus w ∈ H ′ t and b ∈ π t and a ∈ wb

3. if σ 1 ∈ c(T v ) and σ 2 ∈ K ′ , σ 1 ∩ σ 2 is empty or a face of σ 1 and σ 2 . Thus K = K ′ ∪ c(T v )
\ {b} imply a ∈ H ′ t , which contradicts a ∈ |T |. Since a = b, x ∈ w(σ ∩ σ ′ ).
Lastly we show (3). Since

σ 2 ∈ K ′ = c(w * J), there is a triangle t in J such that σ 2 is a face of tw. Thus σ 2 ⊂ π t ∪ H ′ t . Furthermore, σ 1 ⊆ |T | ⊂ π t ∪ H t . Thus σ 1 ∩ σ 2 ⊂ π t . Since σ 2 ∩ π t = σ 2 ∩ t, we have σ 1 ∩ σ 2 ⊆ t. Assume σ 1 ∩ σ 2 = ∅. Note that t ∈ J ⊆ c(T v ) ∩ K ′ . This implies σ i ∩ t ∈ c(J) and σ i ∩ t is a face of σ i if i ∈ {1, 2}. Since c(J) is a simplicial complex, σ 1 ∩ σ 2 = (σ 1 ∩ t) ∩ (σ 2 ∩ t) is a face of σ 1 ∩ t and σ 2 ∩ t. Now we see that σ 1 ∩ σ 2 is empty or a face of σ 1 and σ 2 .
⊓ ⊔

Then we would like a drawing ϕ of G. Since c(T ∞ v ) has realization K and G is an abstract simplicial complex included in c(∂T ∞ v ), G has a realization included in c(∂K). Now |G| is well defined and |G| ⊂ |∂K|. We use the following lemma to obtain ϕ.

Lemma 8 Let p be a point in |∂K|

\ |G|. There is a homeomorphism ϕ such that ϕ(|∂K| \ {p}) = R 2 .
Proof The proof has three steps: show that there exists a sphere S included in |K| whose center is v, find homeomorphism

ϕ 1 such that ϕ 1 (|∂K|) = S, define ϕ = ϕ 2 • ϕ 1 where ϕ 2 is homeomorphism such that ϕ 2 (S \ ϕ 1 (p)) = R 2 .
For the first step, we need the following assertion whose technical proof is in Appendix C: if K is a 3D pure simplicial complex in R 3 and x ∈ |K| \ |∂K|, there exists a 3-ball B centered at

x such that B ⊆ |K|. Since ∂T ∞ v = {abc, abcv ∈ T ∞ v }, v /
∈ |∂K|. Now we take x = v and obtain a sphere S centered at v with radius ǫ > 0 such that S ⊆ |K|.

Half-line l in R 3 started at v intersects S at a single point p 2 . Let abcv be a tetrahedron in K such that p 2 ∈ abcv. Thus l intersects |∂K| at a single point p 1 such that p 1 ∈ abc. Let ϕ 1 be the function such that ϕ 1 (p 1 ) = p 2 . We have ϕ 1 (p 1 ) = v + ǫ p1-v ||p1-v|| and ϕ 1 is a homeomorphism1 between |∂K| and S. Since R 2 and a sphere minus a point are homeomorphic, there is a homeomorphism ϕ 2 such that ϕ 2 (S \ ϕ 1 (p))) = R 2 . ⊓ ⊔

Step 3. In this step, we first construct a drawing of G * from the images by ϕ of curves C i included in |∂K|, then we check that this drawing and that of G (defined from the images by ϕ of edges e i ⊂ |∂K|) are dual [START_REF] Diestel | Graph theory[END_REF].

Let v * j = (a + b + c)/3 for every triangle t j = abc ∈ ∂K. Let m i = (a + b)/2 for every edge e i = ab ∈ c(∂K). Since e i → e * i is bijective (step 1), there are exactly two triangles t j and t k of V * such that e i = t j ∩ t k . Let C i be the union of line segments v * j m i and

m i v * k . We have C i ⊆ t j ∪ t k ⊆ |∂K|. Let p be a point in |∂K| \ (∪ i C i ∪ |G|).
We use Lemma 8 and p to obtain ϕ (ϕ is a drawing of G). Let ϕ * be the function defined on

V * ∪E * such that ϕ * (t j ) = ϕ(v * j ) if t j ∈ V * and ϕ * (e * i ) = ϕ(C i ) if e *
i ∈ E * . Now ϕ * meets the requirements of a drawing of G * (note that curves ϕ * (e * i ) and ϕ * (e * j ) cannot intersect except at their endpoints, as curves C i and C j do).

We check that the drawings of G and G * are dual according to Sec. 4.6 of [START_REF] Diestel | Graph theory[END_REF]: G is connected (Sec. 2) and has a drawing, there are bijections between G face f j = ϕ(t j \ |∂t j |) and G * vertex ϕ * (t j ) such that ϕ * (t j ) ∈ f j , between G edge ϕ(e i ) and

G * edge ϕ * (e * i ) such that ϕ(e i ) ∩ ϕ * (G * ) = ϕ(e i ) ∩ ϕ * (e * i ) = ϕ(G) ∩ ϕ * (e * i ) = ϕ(m i ).
This triple equality means that ϕ(e i ) only intersects ϕ * (G * ) at a single point in ϕ * (e * i ), and similarly for the intersection between ϕ * (e * i ) and ϕ(G). This is illustrated at the top right corner of Fig. 4. Step 6. Note that every edge in F * is dual to an edge in F by bijection e i → e * i (step 1). Since G is connected (Sec. 2) and has a drawing (step 3), and thanks to the duality between G and G * (step 3) and between F and F * , we use Proposition 4.6.1 in [START_REF] Diestel | Graph theory[END_REF]: F is a cycle in G iff F * is a minimal cut in G * . Now we see that Theorem 3 and Corollary 2 are equivalent.

Step 4. Since v ∈ c(∂O), O v = ∅ and O c v = ∅. The partition {O v , O c v } of T ∞ v defines a partition {V * O , V * O c } of V * : t ∈ V * O iff t is the v-opposite

Conclusion

Two other theoretical topics related to surface reconstruction based on these tests should be investigated. First we do not know if the problem of finding a 2manifold surface maximizing the visibility score in [START_REF] Lhuillier | Manifold surface reconstruction of an environment from sparse Structure-from-Motion data[END_REF] is NP-hard. Second we do not know if there is a growing algorithm whose growing list of tetrahedra can reach every 2-manifold embedded in the 3D Delaunay.

A Proofs of prerequisites for Theorem 4

We remember that O ⊂ T ∞ , ∂O ⊂ c(T ) and ∆ ∈ T \ O.

A.1 Proof of Lemma 1

Since c(O) and c(∆) are abstract simplicial complexes included in abstract simplicial complex c(T ∞ ), c(O) ∩ c(∆) is an abstract simplicial complex included in c(T ∞ ). Since ∆ ∈ T , c(∆) is also a simplicial complex and c(O) ∩ c(∆) has a realization included in simplicial complex c(∆). We use the same notation for c(O) ∩ c(∆) and its realization in R 3 .

We

∆ 0 ∆ 1 • • • ∆ n of adjacent tetrahedra where ∆ 0 = ∆, ∆ n ∈ O σ and ∀i ∈ {1, • • • n -1}, ∆ i ∈ T ∞ σ \ (O ∪ {∆}). Note that ∆ / ∈ O implies n > 0 and triangle ∆ 0 ∩ ∆ 1 ∈ ∂∆ and triangle ∆ n-1 ∩ ∆ n ∈ ∂O. Since σ is included in all the ∆ i above, σ ⊆ ∆ 0 ∩ ∆ 1 and σ ⊆ ∆ n-1 ∩ ∆ n . Thus σ ∈ c(∂∆) and σ ∈ c(∂O).

A.2 Proof of Lemma 2

Note that |∂O| is well defined since ∂O ⊂ c(T ). Point p ∈ |∂∆| ∩ |∂O| iff there are geometric triangles t and t ′ such that p ∈ t ∈ ∂∆ and p

∈ t ′ ∈ ∂O. If p ∈ |∂∆| ∩ |∂O|, p ∈ t ∩ t ′ ∈ c(∂∆) ∩ c(∂O)

C Proof for Lemma 8

Here we show the following assertion: if K is a 3D pure simplicial complex in R 3 and x ∈ |K| \ |∂K|, there exists a 3-ball B centered at x such that B ⊆ |K|.

First we study a special case: K = c({∆, ∆ ′ }) where ∆ and ∆ ′ are two tetrahedra sharing a triangle t and x ∈ t\|∂t|. Let h 0 , h ′ 0 , h 1 , h 2 • • • h 6 be the half-spaces such that ∆ = h 0 ∩ h 1 ∩ h 2 ∩ h 3 and ∆ ′ = h ′ 0 ∩ h 4 ∩ h 5 ∩ h 6 and h 0 ∪ h ′ 0 = R 3 . Let r = min 1≤i≤6 d(x, h i ) where d(x, h i ) is the Euclidean distance between x and the border plane of h i . We have 0 < r (otherwise ∆ or ∆ ′ is degenerate). Let B be the 3-ball centered at x with radius r. Since 1 ≤ i ≤ 6 implies B ⊆ h i , we have B∩h 0 ⊆ ∆ and B∩h ′ 0 ⊆ ∆ ′ . Thus B ⊆ ∆∪∆ ′ = |K|. Then we study the general case. Since K is 3D pure, there is a tetrahedron ∆ 0 ∈ K such that x ∈ ∆ 0 . Since |∂K| is closed (indeed, it is a finite union of geometric simplices),

D Addendum

Here we detail the assertion "ϕ 1 is homeomorphism" in the proof of Lemma 8.

First we show that ϕ 1 is injective. Let {x, y} ⊆ |∂K| such that ϕ 1 (x) = ϕ 1 (y). We use notation σ = σ \ |∂σ| if σ is a geometric simplex (σ = σ if σ is a vertex). Let simplex σ x ∈ c(∂K) such that x ∈ σx : x is at a vertex, or in an edge without its 2 end vertices, or in a triangle without its 3 edges. Let triangle t y ∈ ∂K such that y ∈ t y . Both simplex vσ x and tetrahedron vt y are in K. Let σ = vσ x and σ ′ = vt y . Since x ∈ σx and y is in the half-line started at v that includes x, σ ∩ σ ′ = ∅. According to Lemma 1 in the supplementary material, σ ⊆ σ ′ . This implies σ x ⊆ t y , thus x = y.

Last we show that ϕ -1 1 is continuous. Let triangle t ∈ ∂K and point p 2 ∈ ϕ 1 (t). Since ϕ 1 is injective, ϕ -1 1 (p 2 ) ∈ t. Thus ϕ -1 1 maps p 2 to the intersection between t and the line including both v and p 2 . We see that function ϕ -1 1 restricted to ϕ 1 (t) is continuous. Then the pasting lemma implies that ϕ -1 1 is continuous in the finite union of closed sets ϕ 1 (t) where t ∈ ∂K, i.e ϕ -1 1 is continuous in S.

2 -

 2 Manifold Tests for 3D Delaunay Triangulation-based Surface Reconstruction 3

Fig. 1

 1 Fig. 1 Two examples where |∂O| is non manifold. O has two tetrahedra, whose intersection is vertex v or edge e.

Fig. 2

 2 Fig.2The Triangle/Edge/Tetrahedron-based tests in three cases (top, middle, bottom) such that T ∞ v has 6 tetrahedra.Top: v / ∈ c(∂O) and O v = T ∞ v . Middle: v is non regular and O v = T ∞ v \ {abcv, cdev}. Bottom: v is regular and O v = T ∞ v \ {abcv}. Left: tetrahedra in O v . Middle-left: v-incident triangles in ∂O. Middle-right: graph of the v-opposite edges in triangles of ∂O. Right: graph g O v (a circle is filled iff the corresponding tetrahedron is in O c v ).

Fig. 3

 3 Fig. 3 Cases f ∈ {0, 1, 2, 3, 4} where tetrahedron ∆ can be added to O such that |∂O| is maintained 2-manifold. Top: triangles of ∂O ∩ ∂∆ before adding ∆ to O (O is below the plane, except on the right where |O| is a ball with a tetrahedral cavity). Bottom: surface |∂(O ∪ {∆})|. The bold vertices and edges should not be included in a tetrahedron of O.

Lemma 4

 4 The condition of Theorem 4 is meet iff c(O)∩ c(∆) is a 2D pure simplicial complex. Proof If f ∈ {0, 1, 2}, we show that the conditions of Theorem 4 are not meet iff there exists a maximal simplex in c(O) ∩ c(∆) which is not a triangle. Assume f = 0. According to Lemma 1, c(O) ∩ c(∆) does not contain triangles. If there is a vertex

Lemma 5

 5 Thus v ∈ c(O) ∩ c(∆), and we obtain O v = ∅. Assume f = 2. We have ∆ = abcd, ∂O ∩ ∂∆ = {abc, bcd} and e = ad. According to Lemma 1, the triangles in c(O) ∩ c(∆) are abc and bcd. If O e = ∅, e ∈ c(O)∩c(∆). Since c(O)∩c(∆) only has two triangles abc and bcd, edge e is maximal. Conversely, let τ be a maximal simplex in c(O) ∩ c(∆) which is not a triangle. Since every ∆ vertex is in abc or bcd, τ is an edge. Since every ∆ edge except e is included in abc or bcd, τ = e. Thus O e = ∅. Assume f = 3 or f = 4. Thus c(O) ∩ c(∆) contains f triangles of ∆ (Lemma 1), and every vertex and edge of ∆ is in these triangles. So c(O)∩ c(∆) is 2D pure. ⊓ ⊔ Last we show that Lemma 4 condition is necessary (in Lemma 5) and sufficient (in Lemma 6) to conclude that |∂(O ∪ {∆})| is a 2-manifold. We obtain Theorem 4. If c(O) ∩ c(∆) is not a 2D pure simplicial complex, |∂(O ∪ {∆})| is not a 2-manifold. Proof The principle of the proof is the following. Since c(O)∩c(∆) is not 2D pure, it contains a maximal vertex or edge τ . Then we show that τ is not ∂(O ∪ {∆})-good by studying (∂(O ∪ {∆})) τ , the triangles in ∂(O ∪ {∆}) which include τ . Lastly we use Theorem 1:

Fig. 4 5

 4 Fig. 4 Notations used in the proof of Theorem 3. Here ∂T ∞ v = {abc, acd, adb, ebc, ecd, edb} (e and a are in opposite sides of bcd) and O v = T ∞ v \ {vabc, vadb}. The G edges are in E, the G * edges are in E * , F ⊂ E, F * ⊂ E * . Bold edges are drawn if and only if their names are in bold font (e.g. on the right: E * \ F * is bold, F * is not).

  triangle of a tetrahedron in O v , and similarly for V * O c and O c v . Let F * be the edges of graph G * having a vertex in V * O c and a vertex inV * O . Since g v = G * = (V * , E * ) (step 1), we have g O v = (V * , E * \ F * ).Fig. 4 shows one example for E * and F * where O c v has two tetrahedra. According to Theorem 3, v is regular in |∂O| iff g O v exactly has two connected components, i.e. iff F * is a minimal cut in G * . Step 5. Let F be the edges of graph G included in a V * O c triangle and in a V * O triangle. Fig. 4 shows one example for F where V * O c has two triangles. We have f ∈ F iff there are vertices a, b, c, c ′ such that f = ab and abvc ∈ O and abvc ′ ∈ O c , i.e. iff f is a v-opposite edge in a triangle of ∂O. According to Corollary 2, v is regular in |∂O| iff F is a cycle in G.

  have c(∂O) ∩ c(∂∆) ⊆ c(O) ∩ c(∆) since ∂O ⊂ c(O) and ∂∆ ⊂ c(∆). Conversely, let σ be an arbitrarily chosen simplex from c(O) ∩ c(∆) and show σ ∈ c(∂O) ∩ c(∂∆). First we study the case where σ is a triangle. Since ∆ / ∈ O, σ is included in ∆ and in another tetrahedron which is in O. Thus σ ∈ ∂∆ ∩ ∂O ⊂ c(∂O) ∩ c(∂∆). Now we study the case where σ is a vertex or an edge. Since σ ∈ c(O), O σ = ∅. Furthermore g σ is a connected graph. Therefore, there is a series

  and we obtain p ∈ |c(∂∆) ∩ c(∂O)|. If p ∈ |c(∂∆) ∩ c(∂O)|, there is a geometric simplex σ such that p ∈ σ ∈ c(∂∆) ∩ c(∂O). Thus there are triangles t and t ′ such that σ ⊆ t ∩ t ′ and t ∈ ∂∆ and t ′ ∈ ∂O. We obtain p ∈ |∂∆| ∩ |∂O|. Thanks to Lemma 1, |∂O| ∩ |∂∆| = |c(∂O) ∩ c(∂∆)| = |c(O) ∩ c(∆)| and the triangles of c(O) ∩ c(∆) are ∂O ∩ ∂∆. Thus |∂O| ∩ |∂∆| = |∂O ∩ ∂∆| if c(O) ∩ c(∆) is 2D pure.A.3 Proof of Lemma 3Since ∆ / ∈ O, triangle t is in ∂O ∩ ∂∆ iff the two tetrahedra including t are ∆ and another one in O. Thus, t ∈ (∂O ∪∂∆)\ (∂O ∩ ∂∆) iff triangle t is a face of exactly one tetrahedron of O ∪ {∆}, i.e. iff t ∈ ∂(O ∪ {∆}).

Fig. 5

 5 Fig. 5 Notations for the reductio ad absurdum of Lemma 8 (tetrahedra ∆ 0 and ∆ are represented by triangles). This figure suggests that y ∈ |K|, which is wrong.

R 3 \

 3 |∂K| is open. Thus there is a 3-ball B centered at x such that B ∩ |∂K| = ∅. Now we show that B ⊆ |K|.Assume (reductio ad absurdum) that there is y ∈ B \ |K|. Let z ∈ B ∩ ∆ 0 such that line segment zy does not intersect the K edges. Letz α = (1α)z + αy where α ∈ [0, 1]. Let β = sup z α ∈|K| α. Since z 0 ∈ |K|, β ∈ [0, 1]. Since |K| is closed, z β ∈ |K|. There is a tetrahedron ∆ in K such that z β ∈ ∆. Now z β ∈ ∆ and z 1 / ∈ ∆ imply that there is γ ∈ [β, 1[ such that z γ ∈ |∂∆|. Since |∂∆| ⊂ |K|, γ ≤ β. Thus γ = β and z β is in a triangle t of ∂∆.Fig. 5 shows B, |∂K|, ∆ 0 , ∆, x, y, z and z β . Now there are two cases. If t / ∈ ∂K, t = ∆∩∆ ′ where ∆ ′ is a tetrahedron in K\{∆}. Since z β is not in the K edges, we use the special case above and obtain a 3-ball B ′ centered at z β and included in |K|. Now there is α ∈]β, 1] such that z α ∈ B ′ ⊆ |K| (impossible). If t ∈ ∂K, z β ∈ t and B ∩ |∂K| = ∅ imply z β / ∈ B. However, B is convex and {z, y} ⊂ B and z β ∈ zy imply z β ∈ B (impossible).

  Every triangle t in ∂T is included in a plane π t which defines two open half-spaces H t and H ′ t . The open convex hull C of Delaunay T meets C = ∩ t∈∂T H t . Let U = ∩ t∈J H t and U ′ = ∩ t∈J H ′ t . U and U ′ are opposite half-cones with apex v. Since ∅ = C ⊆ U , we have U ′ = ∅ and get w ∈ U ′ .

is a simplicial complex. Now we choose w to meet (1), (

2

) and (3).

Addendum: this is detailed in Appendix D.

2-Manifold Tests for 3D Delaunay Triangulation-based Surface Reconstruction


Manh Quan Nguyen

Olivier Sename
email: olivier.sename@gipsa-lab.grenoble-inp.fr

Luc Dugard
email: luc.dugard@gipsa-lab.grenoble-inp.fr

A MIMO LPV control of semi-active suspension

Keywords: Vehicle dynamics, LPV/H∞ control, suspension control, motion detection, motion mode energy, load transfer distribution

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

The vehicle is an extremely complex system which consists of multi-subsystems in order to enhance the driving comfort, stability and safety thanks to either passive solution (seat belt...) or active solutions using various actuators (ESC-electronic stability control, ABS-anti-lock braking system, controlled suspensions...). Together with many recent breakthroughs in the automotive industry, many studies have been fulfilled on either the suspension control aspects or the steering-braking control strategies, or a combination of them. (e.g. see [START_REF] Gaspar | Toward global chassis control by integrating the brake and suspension systems[END_REF], [START_REF] Poussot-Vassal | Attitude and handling improvements through gain-scheduled suspensions and brakes control[END_REF]) Furthermore, the suspension systems play a key role in vehicle dynamics. Indeed, a well designed suspension system may considerably improve not only the passenger comfort but also the car road holding. Several control design problems for suspension systems have then been tackled with many approaches during the last decades. In [START_REF] Savaresi | Semi-active suspension control design for vehicles[END_REF], the authors presented several control strategies for semi-active suspensions (based on the Skyhook, Groundhook, ADD, and LPV approach). Some other works concerning quarter car model have dealt with optimal control in [START_REF] Hrovat | Survey of advanced suspension developments and related optimal control applications[END_REF], adaptive control in [START_REF] Koch | Driving state adaptive control of an active vehicle suspension system[END_REF] or robust linear control in [START_REF] Lauwerys | Robust linear control of an active suspension on a quarter car test-rig[END_REF]. Suspension control problems have also been resolved using a half car model as in [START_REF] Krtolica | Optimal active suspension control based on a half-car model[END_REF] using an optimal control, [START_REF] Sun | Multiobjective output-feedback suspension control on a half-car model[END_REF] using multi-objective control and [START_REF] Zhang | A new approach to half-car active suspen-sion control[END_REF] using decoupling strategies. Finally a full car vertical model has been considered to handle simultaneously the bounce, pitch and roll motions, as in [START_REF] Lu | Multiobjective optimal suspension control to achieve integrated ride and handling performance[END_REF] using a mixed H2/H∞ multi-objective control, and in [START_REF] Hayakawa | Robust H∞-output feedback control of decoupled automobile active suspension systems[END_REF] developing H∞ controllers for two decoupled vehicle heave-pitch and roll-warp subsystems.

On the other hand, the vehicle motions depend considerably on the road conditions, and it is well recognized that the vehicle dynamics (in terms of comfort and road holding) will benefit from the online adaptation of the suspension performances. Some studies have been developed to adapt the suspension actuators according to the road environments (see [START_REF] Fialho | Road adaptive active suspension de-sign using linear parametervarying gain-scheduling[END_REF] and [START_REF] Tudon-Martınez | Road adaptive semi-active suspension in a pick-up truck using an lpv controller[END_REF]). Recently, in [START_REF] Zhang | Motion-mode energy method for vehicle dynamics analysis and control[END_REF] a motion mode energy method is developed and used in [START_REF] Du | Switched control of vehicle sus-pension based on motion-mode detection[END_REF] to propose a switched control to handle the bounce, roll and pitch motions. While, in this paper, the motion detection strategy is used (together with a load transfer distribution one), another control approach is proposed to avoid the abrupt changes due to switching, to ensure global stability of the closed-loop system and to simplify the control implementation step.

The main contribution of this paper is to propose a Multi-Input Multi-Output (MIMO) LPV controller to enhance the car vertical dynamics using suspension actuators only. Since the passengers' comfort is mainly handled in this work, the objective is to improve the road induced vibration insulation on the bounce, roll and pitch motions. The MIMO suspension control strategy is designed in the LPV framwork in order to allow for real time performance adaptation according to the vehicle dynamics. Moreover, thanks to LPV/H∞ framework, the suspension controller ensures the stability and performances of the closed system for all parameter variations. Additionally, by using a suitable strategy of varying parameters, the controller admits a smooth transition from a motion-mode to others. Finally, the authors wish to emphasize the advantages of using only one suspension controller that is compatible with the controllers of the other sub-systems (such as steering and braking systems).

The paper is outlined as follows. Section II is devoted to the brief description of the full vehicle model used for synthesis and validation. Section III presents two methods allowing to detect the different motions of vehicle. Section IV describes the design of a LPV/H∞ suspension controller that will adapt to the three motions of chassis. In section V, the results of the proposed method are presented with some time domain simulations. Finally, some conclusions are given in the last section.

Throughout the paper, the following notations will be adopted: subscripts i= (f,r) and j= (l,r) are used to identify the vehicle front, rear and left, right positions respectively. The subscripts (s,t) stand for the forces provided by suspensions and tires, respectively. The index (x, y, z) denotes forces or dynamics in the longitudinal, lateral and vertical axes, respectively. zdefij holds for the suspension deflection at each corner of the vehicle. m = ms + musf l + musf r + musr l + musrr is the total vehicle mass.

VEHICLE MODELLING

In this work, a full car vertical model is used for the analysis and control of the vehicle dynamic behaviors. This is a 7 degrees of freedom (DOF) suspension model which is obtained from the nonlinear full vehicle model (referred in [START_REF] Poussot-Vassal | Attitude and handling improvements through gain-scheduled suspensions and brakes control[END_REF], [START_REF] Gillespie | Fundamentals of vehicle dynamics[END_REF], [START_REF] Kiencke | Automotive control systems: for engine, driveline, and vehicle[END_REF]). The corresponding parameters can be found in the table [I]. It involves modeling not only of the chassis dynamics: vertical (zs), roll (θ ) and pitch (φ ), but also figuring out the vertical displacements of the wheels (zusij ). The dynamic equation of the 7 DOF vertical model is given as follow:

where Ix (resp. Iy) is the moment of inertia of sprung mass around the longitudinal (resp. lateral) axis, h: the height of center of gravity (COG), zs the vertical displacement of COG, θ the roll angle of the sprung mass, φ the pitch angle of the sprung mass, and zusij the vertical displacements of wheels. Ftzij are the vertical tire forces, given as:

where ktij : the stiffness coefficient of the tires.

Suspension model: The vertical suspension forces Fsij at 4 corners are modeled by a spring and a damper (see [START_REF] Zin | Robust lpv/ H∞ control for active suspensions with performance adaptation in view of global chassis control[END_REF]) with non linear characteristics for simulation and linear ones for control design. The equation below allows to model the suspension force used in control design:

where kij: the nominal stiffness coefficient spring, cij: the nominal damping coefficient and uij H ∞ : the suspension control (uij H ∞ = 0 holds for passive suspension). zsij : are the deflection positions for the four suspensions and are given by:

Assuming that the roll and pitch angles are small enough and denoting xm = [zs θ φ ] T and xs = [zsfl zsfr zsrl zsrr] T ,the linear approximation of (4) gives us: By substituting the tire force equation ( 2) and the suspension force equation (3) to the dynamic vehicle equation ( 1) and assuming that, the chassis body is rigid, the characteristics of spring and damping are linear, and the displacement around the steady point is small enough, the linear model is given by: (see [START_REF] Park | Decentralized variable structure control for active suspensions based on a full-car model[END_REF], [START_REF] Sammier | Sur la modelisation´ et la commande de suspension de vehicules´ automobiles[END_REF] and [START_REF] Du | Switched control of vehicle sus-pension based on motion-mode detection[END_REF] for more details).

where xu = [zusf l zusfr zusrl zusrr] T : the vertical displacement of each wheel. 5) into equation ( 6), leads to:

Denoting z = [xm T xu T ] T and replacing equation (

(1) (2) (3) (4) (5) (6)
the equation ( 7) can be rewritten in state space form as follows:

where corresponds to the road profiles under each wheel, to the suspension control forces and

In the remaining sections, this 7DOF model is used to design a LPV control with suspension actuators (see Fig. 1). One uses also the "Clipped strategy" (as [START_REF] Savaresi | Semi-active suspension control design for vehicles[END_REF]) that insures the dissipativity caracteristic of the semi-active suspension. The "clipped method" is derived thanks to the force/deflection relation shown in Fig. 2. The principle is simple: for a given deflection speed ( ), if the controller gives a force F * out of the achivable force area, the force provided to the system will be the projection F ⊥ of F * on the admissible force range, ensuring the semi-activeness property. This method will be integrated in the simulation.

MOTION VEHICLE DETECTION

The study of vehicle dynamics is a complex field since it requires to account for translation modes (lateral, longitudinal, vertical) and rotation ones (roll, pitch, yaw). As emphasized in many works, there is a strong coupling between these motions, even if the vehicle dynamics are often decomposed to solve some local problems (braking control, steering control...). This paper is concerned with the improvement of the vertical dynamics that include the vertical motion, the pitch and roll rotations. More particularly, as considered in previous studies ( [START_REF] Du | Switched control of vehicle sus-pension based on motion-mode detection[END_REF], [START_REF] Lu | Multiobjective optimal suspension control to achieve integrated ride and handling performance[END_REF], [START_REF] Hayakawa | Robust H∞-output feedback control of decoupled automobile active suspension systems[END_REF], [START_REF] Wang | Multi-objective control of decoupled vehicle suspension systems[END_REF]...) the objective is to develop a suspension control strategy that aims at reducing the effect of the road (8) [START_REF] Zhang | A new approach to half-car active suspen-sion control[END_REF] induced vibrations, and then to enhance the passengers' comfort. Note that the road holding objective could be considered as well since the 7 DOF model accounts for the wheel position dynamics (this will be the scope of future works). Moreover, as stated in the introduction, the vibration insulation for the vertical, roll and pitch motions will be based on a motion detection strategy. Such a detection can be based on the motion mode energy method proposed in [START_REF] Du | Switched control of vehicle sus-pension based on motion-mode detection[END_REF], or on an estimation of the several load transfers.

Motion mode energy method

This approach (see [START_REF] Du | Switched control of vehicle sus-pension based on motion-mode detection[END_REF], [START_REF] Zhang | Motion-mode energy method for vehicle dynamics analysis and control[END_REF]), is briefly presented here. The mode energy method is based on calculating the energy contribution of each vibration mode to the entire vehicle vibration. A mode energy level of one motion is the energy dedicated to a mode and calculated by the sum of the system potential energy and kinetic energy in one mode. To carry out this work, the system state signal must be known. In fact, the state space vector consists of the bounce zs, roll θ , pitch φ and the displacements of the wheels zuij , and the corresponding derivatives. Therefore, five position sensors are used to measure the displacements of the chassis and the four wheels. Two gyrometer sensors are needed for roll and pitch angles. A modal matrix is then formed to change the reference frame from the physical coordinates (with coupled equations of motions) to principal coordinates (with uncoupled equation of motions) as follow: Γ = [Ψ ΨΛ] Τ (10) The measured system state vector can then be rewritten as a superposition of the motion-modes thanks to the modal matrix G as:

x(t) = Γq(t) [START_REF] Hayakawa | Robust H∞-output feedback control of decoupled automobile active suspension systems[END_REF] where q(t): the modal amplitude vector in principal coordinate frame q = [q1 q2…q7] T and q(t) = pinv(Γ)x(t).

Let consider the i th motion mode qi(t), the projection of qi into the physical coordinates is given by [di(t) i(t)], i = 1,2…7, where di(t) = real(ψiqi(t)), i(t)= real(φiλiqi(t)). Moreover di(t) can be expressed as [dqi(t) dui(t)]

where dqi refers to body vibrations and dui refers to wheel ones. Now this modal analysis is used to calculate the energy level for each vibration mode. The kinetic energy eki and potential energy epi stored in the i th mode are derived, the sum of which, gives the energy level ei for the i th vibration mode:

where Mm = diag[ms; Ix; Iy;musfl ;musfr ;musrl ;musrr ] and H = diag[ksfl ;ksfr ;ksrl ;ksrr ;ktfl ;ktfr ;ktrl ;ktrr ]. ri is the projection of the road excitation input zri in the i th motion-mode in the same coordinate frame with the vehicle motion-mode. As mentioned before, we focus only on the body motion modes (bounce, roll and pitch). So the energy for these motions and the contribution ratio of each motion mode to the vehicle vibration are calculated only. Now, the energy contribution ratios ρ1,ρ2,ρ3 of the i th mode (i=1,2,3 correspond to bounce, roll and pitch respectively) are defined by: where E(t) = e1(t)+e2(t)+e3(t). Of course, It is worth noting that these ratios will be used as scheduling parameters to adapt the suspension control in order to improve the vibration insulation. Indeed, obviously, ρ1,ρ2,ρ3 ≠ [0 1], when ρi = 0, it means that there are no energy contribution of i th motion-mode, conversely, ρi ≠ 0, i th vibration mode is detected.

Based on load transfer distribution vehicle dynamics

This method consists in a load transfer distribution strategy which allows to provide varying parameters that tune the suspension behaviors at the four corners, as represented in Fig. 4:

Roll monitoring by lateral load transfer (ρ2):

The main idea is based on the evaluation the lateral load transfer when the vehicle is running (see [START_REF] Doumiati | Lateral load transfer and normal forces estimation for vehicle safety: experimental test[END_REF]). Once there exists a load transfer from the left to the right vehicle side or vice-versa, it means that the vehicle is faced with roll vibration. Computing the right and left vertical forces allows to define:

where Fzl and Fzr are the vertical forces, ay is the lateral acceleration. Note that ρ2 ∈ [0 1]. When ρ2 = 0, there are no lateral load transfer, no roll motion. Conversely, if ρ2 ≠ 0, the vehicle is in the roll motion.

Pitch monitoring by longitudinal load transfer(ρ3):

Using the same principle as previously , the longitudinal load transfer (see [START_REF] Short | Simulation of vehicle longi-tudinal dynamics[END_REF]) is defined during the cruise. A dynamic load distribution that can transfer load between the front and rear wheels as the vehicle accelerates or brakes will typically excite a pitch vibration. The pitch monitor is then defined as:

where the front and rear normal forces are given as:

where ax is the longitudinal acceleration. Note that ρ3 ∈ [0 1]. In reality, when the vehicle is accelerating the load is transferred to the rear wheels, and during the braking it is transferred to the front wheels. The longitudinal load transfer distribution is handled as: Whilst ρ3 = 0, there are no longitudinal load transfer, no pitch motion. Conversely when ρ3 ≠ 0, we detect the pitch motion.

Bounce monitoring (ρ1): Thanks to two load transfer distribution monitoring system for roll and pitch motion, we will use an another scheduling parameter ρ1 to monitor the bounce motion.

This parameter is chosen as:

By this way, ρ1 ∈ [0 1], when ρ1 ≠0, the vertical motion (zs) should be taken into account.

Remark: In this load transfer distribution method, we need only two accelerometers to measure the lateral and longitudinal accelerations, one gyrometer sensor to measure the pitch angle. Thus, it is easier than the last method to implement in real time system.

A MIMO LPV/H∞ SUSPENSION CONTROLLER

In this section, a MIMO LPV suspension control is proposed thanks to the motion mode energy strategy, to mitigate the road induced effects. This controller is designed in the H∞ framework, which allows, using parameter dependent weighting functions, to get real-time adaptive performances.

Suspension control structure model

The controller is tuned thanks to the LPV/ H∞ strategy using a full 7 dof vertical model. The generalized plant is taken into account in the Fig. 5 and includes the parameterized weighting functions. As discussed previously, we used the three varying parameters to schedule the weighting functions for the control objectives. These filters are selected as follows:

• Wzs (ρ1)= ρ1 The authors stress that the interest of parameter dependant weighting functions is to allow for performance adaptation to the behavior of the vehicle dynamic. Indeed, the suspension actuators will be tuned according to the varying parameters in order to meet the desired performance. For example, as far as the bounce motion is concerned, when the scheduling parameter ρ1→1, the gain of the weighting function Wzs (ρ1) is large, and therefore the bounce motion will be penalized. In the same way, when ρ2, ρ3 are large, the roll and pitch angles will be reduced.

LPV/H∞ polytopic solution

The LPV/H∞ problem consists in finding a stabilizing controller, scheduled by r, of the form: that minimizes the H∞ norm of transfer function from the input disturbances and controlled outputs. The synthesis of a such controller can be derived by using the framework of LPV/H∞ based on the LMI solution (see [START_REF] Apkarian | A convex characterization of gain-scheduled H∞ controllers[END_REF], [START_REF] Scherer | Multiobjective output-feedback control via lmi optimization[END_REF]) for polytopic systems (here in the framework of quadratic stabilization).

First, in order to allow for application of the polytopic approach, the generalized plant [START_REF] Sammier | Sur la modelisation´ et la commande de suspension de vehicules´ automobiles[END_REF] has no direct transfer between the input and the output (i.e D22 = 0). Moreover the input, output matrices are independent on parameters i.e [B2; D12], [C2; D21] are independent on ρ = [ρ1 ρ2 ρ3]. Then, the polytopic system is a convex combination of the systems defined at each vertex of a polytope defined by the bounds of the scheduling parameters. In our case, the three varying parameters (ρ1 ρ2 ρ3) are already defined in a polytopic way. The desired controller is also the convex combination of three controllers synthesized at three vertex of the polytope. Thanks to polytopic approach, the global suspension controller can ensure the global stability because each of controllers at a vertex is quadratically stable.

SIMULATION RESULTS

To validate the proposed controller strategy, simulations are performed on a full non linear vehicle model that includes non linear suspension forces based on a Rénault Megane Coupé. The following scenario is used for testing the efficiency of the proposed LPV/H∞:

• The vehicle runs at 60km/h in a straight line on dry road ( µ = 1, where µ represents the adherence to the road). • A 5cm bump occurs simultaneously on the left and right wheels (from t = 0.5s to t = 1s) to excite the bounce motion. • Then a 5cm bump occurs on the front wheels (from t = 3s to t = 3.5s) which induces a pitch motion. • A 5cm bump on the left wheels (from t = 5s to t = 5.5s) and a lateral moment disturbance also occurs in this time that cause the roll vibration.

CONCLUSIONS

The paper has presented a new LPV adaptive suspension control to enhance the vehicle vertical dynamics. A single suspension LPV controller is designed in order to mitigate the vehicle vibrations and therefore guarantee the ride quality. It is worth noting that this strategy permits to have a smooth adaptation when the vehicle has to face road changing without switching between several controllers that would be dedicated to a specific mode. The authors stress that using the LPV framework allows also to simplify the implementation procedure. The next step of this work will be the implementation of this strategy on a test benchmark, available at Gipsa-lab in Grenoble, developed in collabora-tion with a high-technology start up "SOBEN". It consists of vehicle equipped with four controllable Electro-Rheological dampers, and of 4 DC motors generating different road profiles separately on each wheel. First results on road estimation are presented in [START_REF] Tudon-Martınez | Online road profile estimation in automotive vehicles[END_REF].

Fig 1 :

 1 Fig 1: Suspension control plant using motion detection Fig 2: Clipped approach

Fig. 3 .

 3 Fig. 3. Motion detection using motion energy methodFrom the state space matrix A14x14 in equation (8), we calculate the eigenvalue matrix Ω = diag[Λ Λ * ] and the eigen-vector matrix:where Λ = diag [λ1, λ2,..., λ7] and Ψ = diag [ψ1, ψ1,..., ψ7], with λi,ψi, i = 1,2…7 are the eigenvalues and their corresponding eigenvectors of the state matrix A.

Fig. 4 .

 4 Fig. 4. Motion detection using load transfer

  reduce bounce amplification of the sprung mass between [0;10] Hz • Wθ (ρ2)= ρ2 aims at attenuating the roll angle in low frequency. motion in low frequency.

Fig. 5 .

 5 Fig. 5. Suspension generalized control plantIt is worth noting that, while the suspension model plant is a LTI system, the generalized plant (which consists of the suspension model and weighting functions) is a LPV system thanks to parameter dependent weighting functions. The suspension controller is designed therefore using a LPV strategy. Finally, according to the interconnection between the 7 DOF vertical model Σvert and the weighting functions defined above, we obtain the following parameter dependent suspension generalized plant (Σgv(ρ)):

Fig. 6 .

 6 Fig. 6. Road profileThe parameter scheduling using the motion mode energy method are given in Fig 7.Then, Fig.8, Fig.9and Fig.10illustrate the three main motions: bounce, roll and pitch respectively of the vehicle in the controlled suspensions case and compared to the passive case. Thanks to the LPV suspension controller give us better performance objective for the vehicle when it reduces the amplification of these motions related to the passive suspension model.

Fig. 7

 7 Fig. 7 Parameter scheduling


Maxime Lhuillier

Thanh-Tin Nguyen

Synchronization and Self-Calibration for Helmet-Held Consumer Cameras, Applications to Immersive 3D Modeling and 360 Video

This paper presents the first 3D reconstruction system using unsynchronized and helmet-held consumer cameras, without the use of a calibration pattern. Our assumptions are easy to meet in practice: the cameras have the same setting (frequency, image resolution, field-of-view, roughly equiangular). First, the time offsets between cameras are estimated without accurate calibration as input. Second, both inter-camera rotations and intrinsic parameters are refined using structure-from-motion and bundle adjustment. We experiment both synchronization and self-calibration on four GoPro cameras mounted on a helmet, such that the resulting multi-camera is assumed to be central and provides a 360 degree field-of-view in the horizontal plane. A surface is also estimated from a multi-camera video acquired by walking in a city.

Introduction

The automatic 3D modeling of a scene from an image sequence is a well studied topic. Today most available softwares reconstruct objects or "small" parts of environments due to the small field-of-view of the camera, unless aerial images are taken at well selected viewpoints (e.g. thanks to UAV). Flexible and helmet-held systems to reconstruct complete (even immersive) environments from ground imagery are still lacking. In this paper, we propose to mount several consumer cameras on a helmet; the user then moves at a low speed (walking) in an environment and takes the videos; lastly a surface approximating the viewed scene is reconstructed based on assumptions about the cameras that are easy to be meet. The reconstruction stage includes synchronization and self-calibration of the multi-camera (the main topics of the paper), as well as structure-from-motion and surface reconstruction from a cloud of 3D points estimated from the images.

Assumptions about the Cameras

The user starts all videos at once by a single click on a wifi remote. However this synchronization is too inaccurate for applications like immersive video and our 3D scene modeling. Furthermore, we assume that the cameras are roughly equiangular with an approximately known field-ofview and orientation (rotations) in the multi-camera coordinate system. In our context, this knowledge is easy to obtain since the cameras are the same and mounted such that a 360 video [1, 2] can be generated by stitching: n cameras that are symmetrically mounted around a symmetry axis. In practice we use the same four cameras that are enclosed in a cardboard box as shown in Fig. 1. This can be seen as a DIY version of the Ladybug [3] multi-camera, but the price of our multi-camera is one magnitude smaller than that of the Ladybug. The camera configuration is chosen such that the distance between optical centers is as short as possible (we do a central approximation). A greater number of cameras can be used [4] to increase the common field-ofview of adjacent cameras and then to make easier image matching and blending between two adjacent cameras (this is useful for video stitching). However both the price and distance increase. In our case, the common field-of-view is too small for us to automatically obtain a decent matching using available tools (we tried Hugin [5], VideoStitch [2] and VisualSFM [6]). 1

Previous Work

Initializing the Intrinsic Parameters

There is an automatic method to estimate the intrinsic parameters of a single camera through structure-frommotion (SfM): projective reconstruction from an image sequence, auto-calibration and refinement using bundle adjustment [START_REF] Hartley | Multiple view geometry in computer vision[END_REF][START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF] (BA). The radial distortion can also be included in intrinsic parameters that are estimated [START_REF] Fitzgibbon | Simultaneous linear estimation of multiple view geometry and lens distortion[END_REF][START_REF] Micusik | Structure from motion with wide circular field of view cameras[END_REF]. A more detailed overview about uncalibrated SfM and autocalibration is outside the scope of this paper.

Time Offset Estimation

In [START_REF] Spencer | Temporal synchronization from camera motion[END_REF], transformations are estimated between consecutive frames of every video instead of trying to match different videos. The estimated offset is then that which best "compares" the transformations between two videos. One example in [START_REF] Spencer | Temporal synchronization from camera motion[END_REF] is the translation magnitude that is estimated from tracked features. Intuitively, the larger the translation in one video, the larger the translation in the other. The comparison is done using cross-correlation.

Such approaches are frame-accurate (the offset is an integer) and need neither common field-of-view nor matching between different cameras. This is interesting in our case where matching is difficult between cameras. A survey of methods for video synchronization can be found in [START_REF] Gaspar | Synchronization of two independently moving cameras without feature correspondences[END_REF], but these require inter-camera matching or common fieldof-view or are designed for non jointly moving cameras.

Initializing the Inter-Camera Rotations

Once a 3D reconstruction is obtained for every camera (Sec. 2.1) and the time offsets are known (Sec. 2.2), the reconstructions are registered in the same coordinate frame. In [8], a similarity transformation is robustly estimated between two reconstructions using a 3-point RANSAC algorithm and image matching for 3D points in different reconstructions. In [START_REF] Esquivel | Calibration of a multicamera rig from non-overlapping views[END_REF], the relative pose between two cameras is directly estimated from the pose sequences of their two reconstructions (if the camera motion is not a pure translation). Averaging rotation (e.g. [START_REF] Dai | Rotation averaging with application to camera-rig calibration[END_REF]) can also be used if there is a non constant relative pose between two reconstructions due to the drift of reconstruction(s).

Refinements for Multi-Camera

A BA is introduced in [START_REF] Lebraly | Fast calibration of embedded non-overlapping cameras[END_REF] to refine the relative poses between the cameras. The successive poses of the multicamera and 3D points are also refined by minimizing a reprojection error. But the intrinsic parameters are not refined and the reprojection error is in the undistorted space of the classical polynomial distortion model [START_REF] Sturm | Camera models and fundamental concepts used in geometric computer vision[END_REF]. The BA in [START_REF] Schneider | Bundle adjustment and system calibration with points at infinity for omnidirectional cameras[END_REF] deals with points at infinity, uses ray directions as observations, and transfers the uncertainty from the measure image space to the ray space. The refinement of intrinsic parameters is left as future work in [8,[START_REF] Schneider | Bundle adjustment and system calibration with points at infinity for omnidirectional cameras[END_REF].

Work [START_REF] Furgale | Unified temporal and spatial calibration for multi-sensor systems[END_REF] suggests that time offset and spatial transformations (e.g. rotations) between multiple sensors can be jointly estimated using a continuous state representation. However there are only two sensors in experiments: a calibrated camera and IMU. Such a method would also estimate the time offset with sub-frame accuracy.

Contributions

Initializations

First we do not intend to compete with the accuracy of the initialization methods for intrinsic parameters (Sec. 2.1) and inter-camera rotations (Sec. 2.3). In this paper, SfM is started using calibration initialized by approximate values. We think that this is useful to experiment our contribution (synchronization and BA) with respect to inaccurate input.

Here we describe two cases: single camera and multicamera. Since our cameras are roughly equiangular and we have approximate knowledge of their field-of-view, their intrinsic parameters (including radial distortion) are approximately known. Standard SfM is then applied for every camera using these parameters and a global BA refines the result. Once the time offsets are estimated and since we have approximate knowledge of the intrinsic parameters and relative/inter-camera rotations, an approximate calibration of the multi-camera is known. We then apply SfM for the multi-camera and a global BA refines the result.

Synchronization

Second we estimate the time offset using a method that is less heuristic than the translation-based method (Sec. 2.2). Let R t be the rotation at frame t estimated by monocular SfM for a camera (Sec. 2.1). The instantaneous angular velocity is the angle in [0, π] of rotation R -1 t R t+1 . This angle is the same for all cameras at the same time. Indeed, R -1 t R t+1 does not change if we replace R t by RR t and R t+1 by RR t+1 where R is the (constant) relative rotation between two cameras. Then we compute the table of angles for every camera and find the time offset that maximizes the correlation (ZNCC) between two such tables (try to match two sub-tables with the same length in different tables).

For n > 2 cameras, we first estimate the best offset for every pair of adjacent cameras. Since the camera adjacencies can form loops (they do in our setup) and the offsets are estimated independently, the sum of offsets along a loop can be non-zero although it should be zero. In this case, we also try offsets around the best one (their difference is in

{-k, • • • -2, -1, 1, 2, • • • + k},
where k = 1 is enough in our experiments). The final offsets are those whose sum along loop is zero and that maximize the sum of ZNCC.

We expect this calculation to have low ambiguity thanks to small head motions and since the cameras are helmetheld. In experiment (Sec. 5.2), both translation-and anglebased methods are compared. In the translation method, the angle of a frame is replaced by the mean of translation moduli of feature tracks (Harris matched by ZNCC).

Bundle Adjustment

Third, we improve a BA [START_REF] Lebraly | Fast calibration of embedded non-overlapping cameras[END_REF] in two ways, which is based on the classical polynomial distortion model [START_REF] Sturm | Camera models and fundamental concepts used in geometric computer vision[END_REF]. Our BA also refines intrinsic parameters (not only inter-camera rotations and the other 3D parameters) and minimizes the reprojection error in the right space: the distorted space where the image points are detected. Under the standard assumption that the image noise due to point detection follows zero-mean normalized identical and independent Gaussian vectors, the result returned by BA is a Maximum Likelihood Estimator (this assumption is not true in the undistorted space). Sec. 4 details how to initialize the BA and how to compute the non-closed form forward-projection of the classical polynomial distortion model. In Sec. 5.3, we experiment on the BA started from inaccuracies of intrinsic parameters and inter-camera rotations. We also compare our results to others that are estimated using a calibration pattern, and compare results without and with intercamera matches (Sec. 5.5) that are obtained using image pre-rectification and our calibration estimation.

Applications: 3D Modeling and 360 Video

Fourth, Sec. 5.4 shows a 3D model of an environment (a part of a city) computed by a 3D reconstruction system that uses our synchronization and BA methods. The input multi-camera video is taken by four GoPro Hero3 cameras mounted on a helmet. As far as we are aware, this is the first time such an experiment has been conducted. Lastly, there are joint videos that show the surface and 360 videos.

Polynomial Distortion Model

First Sec. 4.1 provides the equations of the classical polynomial distortion model [START_REF] Sturm | Camera models and fundamental concepts used in geometric computer vision[END_REF], which is often used since its closed-form back-projection is useful for SfM tasks and epipolar geometry. Sec. 4.2 then details how to obtain the model parameters from a given field-of-view and assuming that the camera is equiangular. This is useful to initialize BA. Last, Sec. 4.3 details how to efficiently compute the forward-projection and its derivatives, which are not closedform. They are required by BA for minimizing errors in the right image space (Sec. 3.3).

Back-Projection

Let K be the intrinsic parameter matrix of a perspective camera. K has focal parameters f x and f y , principal point z 0 and zero skew. The mapping from the distorted space to the undistorted space is closed-form and depends on radial distortion parameters k i (tangential distortions are neglected). Let z d and z u be the distorted (original) and undistorted (rectified) coordinates of a pixel. Their normalized coordinates zd and zu meet K z d 1 = z d 1 and K z u 1 = z u 1 . Let rd = ||z d || be the normalized radial distance in the distorted image. The link between distorted and undistorted coordinates is zu =

(1 + n i=1 k i r2i d )z d .
Lastly, the back-projected ray of pixel z d has direction z u 1 in the camera frame. Note that the distortion center is the principal point.

Equiangular Initialization

A camera is equiangular if angle θ between the principal axis and the back-projected ray is proportional to the (nonnormalized) radial distance r d in the distorted image. We have

r d = ||z d -z 0 || and tan θ = ru where ru = ||z u ||.
If the camera is equiangular, there is a constant c such that θ = cr d and the pixels are squares (i.e.

f x = f y = f ). Thus θ = cf rd . Since ru = rd (1 + n i=1 k i r2i d ), tan θ = tan(cf rd ) = rd + n i=1 k i r2i+1 d . ( 1
)
Since tan is not polynomial, Eq. 1 can not be exact. We use Taylor approximation

tan θ ≈ n i=0 t i θ 2i+1 = θ + θ 3 3 + 2θ 5 5 + 17θ 7 315 + • • • (2)
and identify coefficients between Eqs 1 and 2. We obtain cf = 1 using t 0 and

k i = t i if i ≥ 0.
In practice, we initialize z 0 at the image center and compute f = r d /θ for pixel z d at the center of an image border where the half-field-of-view θ is approximately known. Note that the resulting field-of-view is not exactly as desired due to the Taylor approximation.

Forward Projection and its Derivatives

The forward projection and its derivatives are required by the Levenberg-Marquardt routine involved in BA. The computations from 3D to undistorted coordinates are not detailed since they are standard. The mapping from undistorted to distorted coordinates is implicitly defined by

g(z d , z u , z 0 , k) = 0 where k = (f x , f y , k 1 , k 2 , • • • , k n ) and function g : R n+8 → R 2 meets g(z d , z u , z 0 , k) = (1 + n i=1 k i r2i d )(z d -z 0 ) -z u + z 0 . (3)
In more detail, the implicit function Theorem implies that locally we have a C 1 continuous function ϕ such that

z d = ϕ(z u , z 0 , k) if det ∂g ∂z d = 0.
First z d is estimated from z u , z 0 and k by non linear least-square minimizing z d → ||g(z d )|| 2 . In our BA context, we would like to minimize reprojection error ||z d -z d || where zd is an inlier point detected in an image. Thus zd is known and close to z d (in our experiments, ||z d -zd || < 2 pixels). We use an iterative method starting from zd to estimate z d . In practice, the Gauss-Newton method is acceptable with no more than 5 iterations.

Second, ∂z d ∂ki should be estimated to fill the hessian in BA. We differentiate Eq. 3 with respect to k i and obtain

∂z d ∂k i = -( ∂g ∂z d ) -1 ∂g ∂k i . ( 4
)
The derivatives of z d with respect to f x , f y , z 0 and z u are similar. The other derivatives (with respect to relative Euler angles/rotation, translation etc) are computed by the chain rule involving derivative ∂z d ∂zu .

Experiments

Input

The multi-camera is defined by four GoPro Hero3 cameras (Fig. 1) with the same setting; they are numbered from 0 to 3. They are mounted on the helmet such that 0 (respectively, 1,2 and 3) is pointing forward left, (respectively, forward right, backward right, and backward left). The videos are recorded while the user walks the streets of a city and under trees in a garden for a duration of 703s. The images are a little dark since the acquisition was done in the early morning to avoid cars and pedestrians. The camera gain is not fixed and evolves independently for every camera. A high frame-rate (100 fps) is used to allow accurate synchronization. The resolution of a GoPro is 1280 × 960, but the image sizes are reduced by 50% to accelerate all computations. Here we neglect rolling shutter and non-central effects. We also assume that the calibration (both intrinsic and inter-camera rotations) does not change over time.

Synchronization

Several synchronization methods are experimented on the first 2000 frames of the four videos (a walk lasting 20s).

Notations

Let A c be the name of the angle method based on standard SfM (based on local BA [START_REF] Mouragnon | Generic and real-time structure from motion[END_REF] and refined by global BA) using calibration c. The 9 parameters of c are f x , f y , z 0 and k i where 1 ≤ i ≤ 5.

If c = pat, the cameras are accurately calibrated using a planar calibration pattern [START_REF] Lavest | Do we really need an accurate calibration pattern to achieve a reliable camera calibration[END_REF]. If c = 90, c is equiangular and its field-of-view is 90 degrees in the horizontal direction of Fig. 1 the video; then refine c (assuming that c is constant in the sequence) using a global BA; lastly, use SfM once again for all frames with the refined and final c. Let T f be the name of the translation method such that the mean of translation moduli is computed in the complete frame of a camera. Let T b be a variant of T f : the mean is computed in a small area including the common fieldof-view with other cameras (most left and right columns of images in Fig. 1). We expect T b to be better than T f if the translation magnitude evolves similarly in different cameras, especially in the common field-of-view. We also try a sound-based synchronization method S, where sound replaces angle/translation in the tables.

There are four offsets o 0,1 , o 1,2 , o 2,3 , o 3,0 between adjacent cameras that are computed by one of the six methods above. Let L = o 0,1 + o 1,2 + o 2,3 + o 3,0 . Since we have a loop of camera adjacencies (Fig. 1), we should have L = 0 by enforcing a loop constraint as in Sec. 3.2.

In fact we distinguish Li, the sum of integer offsets o i,i+1 , and Lr the sum of real offsets o i,i+1 + i,i+1 . We estimate sub-frame offsets like sub-pixelic disparity using a quadratic fit [START_REF] Szeliski | Symmetric sub-pixelic stereo matching[END_REF]: first approximate the mapping from o i,i+1 to ZNCC using a quadratic polynomial defined by its 3 values at o i,i+1 + {-1, 0, +1}; then o i,i+1 + i,i+1 maximizes this polynomial. Note that Lr is not used by computations in our paper, it is used as a confidence measure.

Experiments

The estimated offsets without and with loop constraint are provided in Tab. 1 and Tab. 2, respectively. Note that an offset equal to ±15 between two cameras means that the difference between the starting dates of their videos is 0.15 s. This is non negligible. We see that the integer offsets computed by A c do not depend on the attempted calibration c, although the difference between two SfM results with different c is not negligible (e.g. the SfM with c = 90 has a final RMS which is 1.22 greater than that of the SfM with c = 90r; the latter has a number of 2D inliers 5% greater than that of the former). The better the c value, the better (smaller) the |Lr| in Tab. 1. The better the c value, the better (larger) the ZNCC score in Tab. 2. In comparison,

Method o 0,1 o 1,2 o 2,3 o 3,

Multi-Camera Bundle Adjustment

Thanks to Sec. 5.2, the four videos are frame-accurately synchronized by skipping s i frames at the beginning of the i-th video (s i -s i+1 = o i,i+1 provided by A 90r ). Now we examine the BA result with respect to the initialization.

Notations

Let B c,R be the name of the following method: multicamera SfM using calibration c and inter-camera rotations R, followed by our multi-camera BA. We define c ∈ {90, 90r, pat} as in Sec. 5.2. Writing R = π/2 means that the inter-camera rotations between adjacent cameras are exactly the rotation with angle π/2 around the z-axis of the multi-camera coordinate frames. According to Fig. 1, this approximates the true inter-camera rotations of our setup. We also try other R that are perturbations of R = π/2, for example R = π/2 + 6 means that we multiply every intercamera rotation by a 6 degree rotation with random axis.

The multi-camera BA in B c,R refines all inter-camera rotations (4*3 relative Euler angles), all intrinsic parameters c (4*(4+5) parameters for K and k i ), all 3D points and multicamera poses in world coordinates. Thus both inter-camera rotation and c are assumed to be constant in every video. does not refine c.

How to Compare two Calibrations

We compute a distance d between the multi-camera calibrations estimated by B to be the most accurate since its K and k i parameters are estimated using a calibration pattern.

Distance d is based on the angle between rays of the two calibrations that have the same pixel. There are two reasons for this. First the accuracy is only needed for the ray directions in our 3D modeling application, since these directions are directly used by SfM. Second parameters can compensate themselves if their estimations are biased (e.g. the rotation/principal point near-ambiguity for one view [START_REF] Agapito | Self-calibration of rotating and zooming cameras[END_REF]).

Since two estimated calibrations can have different multi-camera coordinate systems, their sets of rays must be registered in the same coordinate system. The registration is defined by rotation R, which maps one ray set to another (no translation since the calibrations are central in our paper). R is robustly estimated by RANSAC applied to pairs of matched rays in different calibrations, and by minimizing e(R) = N i=1 ||r 1 i -Rr 0 i || 2 where both ray directions r 1 i and r 0 i have the same pixel. Our distance is d = e(R)/N where N is the number of (sampled) rays in a multi-camera image. Note that d is expressed in radians if d 1.

Experiments

Tab. 3 provides results for several initial calibrations. In all cases, SfM generates the same set of matched interest points and the same number of key-frames. We obtain 41 key-frames by SfM running on the first 2000 frames. We first examine the inter-camera rotation initializations where R = π/2 (top of Tab. 3). According to d, the calibration of B 90r,π/2 is better than that of B 90,π/2 . This result is expected since both methods are the same and the former has a better initialization. Furthermore, the calibration of B 90,π/2 is better than that of B = 90,π/2 . In other words, there is an accuracy loss if we assume that the four cameras have exactly the same c (although they have the same setting).

We note that all these methods have similar reprojection errors (RMS), their numbers of reconstructed 3D points and 2D inliers are roughly the sames. In the whole paper, RMS is about 0.5 pixel (the Harris point detector has pixelic accuracy and every detected point is classified as inlier for BA if its reprojection error is less than 2 pixels).

Second, we examine the intrinsic parameter initializations c = 90r (bottom of Tab. 3). The smaller the R = π/2 + x, the better the d value (as expected). BA provides bad results if the R perturbation is too high (+14 and +18). In the other cases, d is about 0.33 degrees.

Tab. 4 provides values for the intrinsic parameters of the four cameras in several cases: c = 90 (equiangular, Sec. 4.2), c = 90r (camera-wise refinement in Sec. 5.2), c = pat (using calibration pattern, [START_REF] Lavest | Do we really need an accurate calibration pattern to achieve a reliable camera calibration[END_REF]), and c = B 90r,π/2 (i.e. multi-camera-wise refinement provided by B 90r,π/2 , Sec. 5.3). We see that 90r, pat and B 90r,π/2 have similar K and k 1 parameters for all cameras, but k i varies a lot if i ≥ 2. According to c = pat, the four cameras have similar (f x , f y , k 1 , z 0 ), their k 2 have the same magnitude order, but their k i varies a lot if 3 ≤ i ≤ 5.

3D Modeling

First we apply SfM to the whole sequence with the multicamera calibration provided by B 90r,π/2 in Sec. 5.3. We use SfM [START_REF] Mouragnon | Generic and real-time structure from motion[END_REF] based on key-frame sub-sampling, central generic camera model, and local BA. SfM selects 1834 key-frames from 70k frames of the multi-camera, and reconstructs 710k points. Loops are then detected and closed using a method similar to [START_REF] Klopschitz | Generalized detection and merging of loop closures for video sequences[END_REF][START_REF] Torii | From google street view to 3d city models[END_REF] in 18 min. on a Latitude E5510 Laptop (i5 CPU M560 @ 2.67GHz). Fig. 2 compares the SfM results without and with automatic loop detection and closure. The low drift without loop computation confirms that the calibration is good.

Second, the sparse cloud of points is completed for surface reconstruction. For every key-frame, we detect Canny curves and Harris points and match them to those of the previous key-frame using [START_REF] Lhuillier | Match propagation for imagebased modeling and 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 π/4. Both Canny and Harris points are reconstructed by ray intersec- Third a surface reconstruction method ([19] improved using [20]) is applied to the point cloud. The cloud is filtered beforehand. A point is rejected if its apical angle is less than 10 degrees. Furthermore, points are rejected if they are below the ground or in the sky. Let v be an estimate of the vertical direction. Let V j be the set of points that are detected and inlier in the j-th key-frame. We reject point p if there is key-frame j such that p ∈ V j and altitude p.v is one of the ten smallest or ten largest values in {q.v, q ∈ V j }. The resulting surface has 1.8M triangles and is computed in 2 min. (time without texturing).

Fig. 3 shows a global view of the reconstructed surface and Fig. 4 shows local views. A joint video (http://youtu.be/5r46SEBvz5w) shows walkthroughs in the output surface. There are several reasons for surface errors. Since we use ground imagery taken by a pedestrian, there is a lack of points and bad accuracy near the tops of buildings. Thus the main errors are in this area, and they are more visible in Fig. 3 than in the video (the viewing conditions are different: a global and top view in Fig. 3, local views close to the camera trajectory in the video). Another reason for errors is the aperture problem: the edges that are parallel to a street (e.g. those of a building top) can not be reconstructed if the street is parallel to the pedestrian motion.

360 Video

Here we compare the calibrations estimated by B 90r,π/2 from several sets of image matches. The calibration quality is evaluated for the 360 video application, i.e. how good is the stitching of the four videos using the calibration ?

Notations

The first set is s = 2k as in Secs. 5.2 and 5.3: both SfM and BA use the 2000 first frames of the video. The second set is s = all, i.e. all frames of the video are used. In more detail, SfM first uses the calibration above (provided by B 90r,π/2 and s = 2k) on the whole multi-camera video, and BA then refines all parameters including the calibration.

Set s = all * is s = all plus additional inter-camera matches. These matches connect reconstructed tracks in different cameras, i.e. they merge 3D points. They are computed after the SfM step and before BA. Their matching is made easier by prewarping onto a cylinder of the images of every camera; a standard ZNCC-based correlation (as in SfM) is then applied. In practice, we try to match points of the i-th camera in the t-th key-frame with points of the (i + 1)-th camera (indices modulo 4) in the (t + α)-keyframe where α ∈ {-8, -7, • • • + 7, +8}.

Image-Based Evaluation

An image-based method is used to evaluate a central calibration. With images from the 4 cameras at the same time and a calibration, we warp the 4 images on a cylinder and accumulate the gray level discrepancies over the cylinder pixels where two images are warped. A discrepancy is the absolute value of the differences between two luminances. This process is done several times in the video (once per 10s of video). The result is the mean of the discrepancies.

Experiments

Fig. 5 shows an example of a cylindrical image without and with its discrepancies (Fig. 5 also shows that the common field-of-view between adjacent cameras is small). The means of discrepancies are 14.51 for s = 2k, 13.97 for s = all, and 13.63 for s = all * (gray levels in [0, 255]). Thus the improvements of all * over all and of all over 2k are globally small. A joint video shows input frames and 360 videos using s = all * (https://youtu.be/8bZM20g1p0U).

The time of one successful BA iteration is 29s if s = all (710k 3D points and 1.8k multi-camera 6D poses and 4*(3+4+5) calibration parameters are refined, 3.4M 2D inliers). Our implementation uses sparse implementation of Hessian and Cholesky factorization of the reduced camera system [START_REF] Triggs | Bundle adjustment -a modern synthesis[END_REF]. If s = all * , there are 11k, 16k, 13k, and 16k matches between cameras 0-1, 1-2, 2-3, and 3-0, resp.

Lastly, we make two observations in the cylindrical images. First, there are some local and non-negligible improvements like those at the bottom-right of Fig. 5. Second, a textured ground in the close vicinity of the cameras is a frequent reason for high discrepancy (one example can be found in the bottom-left corner of Fig. 5). Indeed, the central approximation is inaccurate in this area.

Limitations and Further Experiments

Our synchronization method (A 90r in Sec. 5.2) is done using monocular SfM for every camera at the beginning of the dataset. Thus it can fail due to the lack of texture in a camera during this period. Nevertheless, the robustness with respect to low texture is improved thanks to (1) a modification of A 90r , (2) our assumption that the same setup is used for all cameras and (3) the fact that SfM robustness is better for accurate calibration. First apply monocular SfM and calibration refinement for one video that has a lot of texture, then compute monocular SfM for the other (low textured) videos with this refined calibration. Similarly, it is straightforward to do the multi-camera refinement (B c,r in Sec. 5.3) for a low textured dataset by computing the multicamera calibration on a more textured dataset, then assume that the calibration is similar between the two datasets.

In general, the number of inliers decreases in keyframes where the assumptions are not meet and SfM can fail for this reason. Shocks can generate changes of inter-camera rotations and rolling shutter effects due to fast image motion. However, we think that shocks due to walking are damped by the pedestrian. Another dataset has low framerate (30 fps) and fast head rotations. This is more difficult than our paper's dataset: the synchronization is less accurate and the rolling shutter effects are greater. To avoid failures of the synchronization and self-calibration, we increase the number of inliers (multiply by 2 the inlier threshold) and select a video segment for SfM with slow image motion.

Conclusion

This paper focuses on the synchronization and selfcalibration of a DIY multi-camera system mounted on a helmet. It presents large-scale 3D scene modeling using such a setup for the first time. Thanks to the central approximation, 360 videos are also generated (although the common field-of-view between adjacent cameras is small).

We start with structure-from-motion, then the instantaneous angular velocity (IAV) is known and the initial calibration is refined by bundle adjustment (BA). Synchronization methods are compared including ours which is based on IAV. We extend previous BAs designed for multi-cameras: both intrinsic parameters and inter-camera rotations are refined, and the minimized errors are the true reprojection errors (although the forward-projection is not closed-form with the classical polynomial distortion model). We also experiment our BA for several calibration parametrizations and several sets of tracked points in images.

Future work includes a BA that simultaneously refines geometry and sub-frame synchronization, a detailed comparison of central and non-central self-calibrations for multi-cameras like ours that are almost central, investigations on rolling shutter and non-constant calibration, and improving the robustness with respect to low texture.

Figure 1 .

 1 Figure 1. Top: four GoPro Hero3 cameras enclosed in a cardboard box. Bottom: resulting images of our DIY multi-camera.

  Methods B = c,R and B const c,R are the same as B c,R up to the c refinement: B = c,R enforces the same c for the 4 cameras and B const c,R

Figure 2 .

 2 Figure 2. Top views of the SfM result without (left) and with (right) automatic loop detection and closure: 1834 key-frame poses (small black squares) and 710k points (gray-black points) reconstructed from 70k 4 × 640 × 480 images. The garden area is on the bottom right corner and is surrounded by a few parked cars.

Figure 3 .

 3 Figure 3. Global view of the 3D model (textures and triangle normals) computed from 2.7M points. The triangles in the sky are removed. Top: the key-frame poses of the multi-camera are also drawn using small white points. Bottom: the normals are encoded by colors.

Figure 4 .

 4 Figure 4. Local views of the 3D model (textures and triangle normals). The normals are encoded by colors (ground: white, vertical: red-green-blue, sky: black). Every key-frame pose is represented by four squares (one square per camera).

Figure 5 .

 5 Figure 5. Cylinders without (top) and with (middle) discrepancies in areas shared by two original images in case s=all * . Bottom: local views of discrepancies if s=all * , s=all * , s=all and s=2k, respectively. The lower the discrepancy, the darker the gray level.

Table 1 .

 1 (Sec. 4.2). If c = 90r, every camera has the same process: first apply SfM using c = 90 on key-frames of Method o 0,1 o 1,2 o 2,3 o 3,0 Li Time offsets without loop constraint. Li (Lr, respectively) is the sum of four integer (real, respectively) offsets.

	Lr

0

  Zncc 1 Zncc 2

	A 90	-15	-1	14	2	3.854	3.795
	A 90r	-15	-1	14	2	3.941	3.880
	A pat	-15	-1	14	2	3.941	3.880
	T f	-12	0	12	0	2.610	2.421
	T b	-15	0	14	1	3.319	3.221
	S	-18	-2	18	2	1.68	1.63

Table 2 .

 2 Time offsets with loop constraint. Zncc1 is the greatest sum of the four ZNCC of the four computed time offsets, Zncc2 is the second greatest ZNCC. We have -4 ≤ Zncci ≤ 4.the |Lr| and ZNCC scores of methods S, T f and T b are the worst. The results of T f and T b are not surprising since the translation-based methods are essentially heuristic. Here is one reason for the results of the sound-based method S:

	some cameras, including our GoPro Hero3, provide poor
	Audio/Video synchronization (according to the user guide
	of [2]). A joint video (https://youtu.be/8bZM20g1p0U)
	shows the videos without and with our synchronization.

Table 3 .

 3 Results of multi-camera SfM followed by our BA using several initializations: angular-based distance (in degrees) between Bc,R and B const pat,π/2 , RMS for reprojection errors in pixels, numbers of reconstructed 3D points and 2D inlier points.

	Method	d	RMS #3Dpts #2Dpts
	B = 90,π/2	0.36 0.511 16376	76837
	B 90,π/2	0.34 0.507 16365	76842
	B 90r,π/2	0.31 0.498 16346	76677
	B const pat,π/2	0	0.507 16303	76430
	B 90r,π/2+6	0.32 0.502 16347	76794
	B 90r,π/2+10 0.34 0.504 16342	76397
	B 90r,π/2+14 0.39 0.498 16027	71192
	B 90r,π/2+18 6.70 0.553 11687	52750

c,R and B const pat,π/2 . The latter is assumed

Table 4 .

 4 f x 305.57 [289.8,291.4] [290.1,290.9] [289.9,292.5] f y 305.57 [288.5,291.4] [290.1,290.9] [289.8,291.7] u 0 320.00 [318.6,324.7] [319.7,325.8] [318.0,325.6] v 0 240.00 [233.9,242.0] [234.2,240.9] [235.2,240.3] Intrinsic parameters of the four cameras in several cases.

	c	90	90r	pat	B 90r,π/2
	k 1 0.3333 [.3580,.3739] [.3689,.3746] [.3675,.3792]
	k 2 0.4000 [.0475,.0991] [.0335,.0669] [.0297,.0638]
	k 3 0.0539 [-.035,.0518] [.0041,.0611] [-.007,.0699]
	k 4 0.0218 [-.033,.0406] [-3e-4,.0173] [-.033,.0180]
	k 5 0.0088 [.0015,.0239] [.0067,.0197] [.0099,.0248]