Luis Rodríguez

Fernando Cuartero

Vicente Pascual

A new approach to the Maximum Clique problem for large random graphs

The search for the maximum clique in a graph is one of the most widely studied problems in graph theory and combinatorial optimization. In the last decades, numerous algorithms and methods to address this problem have been proposed. Most of these algorithms are based on a backtracking scheme to find the maximum clique. In this paper we propose a novel method where the main problem is decomposed in smaller sub-problems which involve finding maximum cliques for certain subgraphs of the given graph. Although these sub-problems are finally solved by using state-of-the-art methods, this decomposition can reduce the complexity of the main problem leading to faster algorithms. Besides, by solving all these sub-problems we obtain additional information about the cliques in the input graph.

Introduction

Given an undirected graph G = (N, E) where N = {1, . . . , n} is the set of nodes and E ⊆ N × N is the set of edges, we can define a clique as a subset U ∈ N where ∀i, j ∈ U : i = j ⇒ (i, j) ∈ E. That is, for any two nodes in U , there is an edge connecting those nodes in E. A clique is also called a complete subgraph.

If U is a clique of G we can say that U is a maximal clique if we can not obtain a larger clique U by adding new nodes from G to U . Similarly, we say that U is a maximum clique of G if U is the largest of all the maximal cliques in G.

The applications of the maximum clique problem are many, ranging from Chemistry to Telecommunications or Artificial Intelligence.

State of The Art

The maximum clique problem is a hard problem since the associated decision problem, which involves determining the existence or not of a k-size clique in a graph has proven to be N P -complete [START_REF] Karp | Reducibility among Combinatorial Problems[END_REF]. As a result of this, we know that, unless P = N P , any algorithm able to solve this problem has, in the worst case, a time complexity which is exponential with the size of the graph.

Several approaches have been proposed to deal with the maximum clique problem. A nice description of the main ones can be found in [START_REF] Bomze | The Maximum Clique Problem[END_REF], and in [START_REF] Wu | A review on algorithms for maximum clique problems[END_REF] we can find an updated and comprehensive review on the algorithms existing in the literature coping with this problem.

Although most of the proposals are quite interesting from a theoretical point of view, the most successful methods in practice are based on a particular algorithm, developed in 1973 by Dutch scientists Joep Kerbosch and Coenraad Bron, [START_REF] Bron | Algorithm 457: Finding all cliques of an undirected graph[END_REF], which relies on a backtracking scheme to explore all the maximal ISSN subm. to DMTCS c by the authors by the author(s) Distributed under a Creative Commons Attribution 4.0 International License cliques in a graph. From this, the maximum clique can be obtained just by returning the largest of these maximal cliques. In [START_REF] Bron | Algorithm 457: Finding all cliques of an undirected graph[END_REF], the algorithmic complexity for the worst case is determined to be in O(3 n/3 ), where n is the number of nodes in the graph. For random graphs, however, no theoretical characterization is so far available, although experimental results show, in general, a performance significantly better than this worst case.

The backtrack Bron-Kerbosch algorithm has been the base of a wide range of methods by using different branch and bound techniques to improve it. In [START_REF] Johnston | Cliques of a graph-variations on the bron-kerbosch algorithm[END_REF], [START_REF] Balas | Finding a maximum clique in an arbitrary graph[END_REF], [START_REF] Carraghan | An exact algorithm for the maximum clique problem[END_REF], [START_REF] Babel | A branch and bound algorithm for the maximum clique problem[END_REF] and [START_REF] David R Wood | An algorithm for finding a maximum clique in a graph[END_REF] among others, we can found a wide set of different algorithms.

In [START_REF] Tomita | An efficient branch-and-bound algorithm for finding a maximum clique[END_REF], an approximated vertex coloring method is used as prune method within this scheme. This coloring algorithm is used to prune those branches that do not improve the largest clique found so far. Furthermore, the coloring algorithm defines a vertex order that is used to first explore those branches (nodes) that seem to be more promising for finding the largest clique. An improved coloring method along with the use of dynamic bounds was proposed in [START_REF] Konc | An improved branch and bound algorithm for the maximum clique problem[END_REF]. In addition, in [START_REF] Tomita | An efficient branch-and-bound algorithm for finding a maximum clique with computational experiments[END_REF] and [START_REF] Tomita | A simple and faster branch-and-bound algorithm for finding a maximum clique[END_REF] several improvements to the proposal in [START_REF] Tomita | An efficient branch-and-bound algorithm for finding a maximum clique[END_REF] can be found. These are based on performing an initial sorting of the graph according to the degree of its nodes along with an enhanced coloring method that tries to reduce the number of colors assigned to the graph nodes.

Recently, in [START_REF] San | An improved bit parallel exact maximum clique algorithm[END_REF] a different graph encoding based on bit string was proposed, which result in a significant speed improvement with respect to the MCS Tomita's algorithm presented in [START_REF] Tomita | A simple and faster branch-and-bound algorithm for finding a maximum clique[END_REF].

The HV algorithm for the maximum clique problem

As described in section 2, the most successful methods to solve the maximum clique problem are based on a backtracking algorithm. In this work, we propose a novel approach where the problem of finding the maximum clique on an input graph is split into several, smaller sub-problems. This strategy is based on a heterodox view of the divide-and-conquer method, where the sets of sub-problems obtained by breaking down the original problem might not be disjoint. Accordingly, some sub-problems may be solved more than once, which results into a less effective reduction process. Thus, the original problem can not be reduced to a polynomial time algorithm, which is consistent with the N P -complete nature of the maximum clique problem.

Formal framework

Let G = (N, E) be the input graph where N = {1, . . . , n} is the sorted set of nodes. Let i, j ∈ N be two nodes where i < j. We can define a subgraph

G ij = {N ij , E ij } where N ij = {k ∈ N / i < k < j} and E ij = {(k 1 , k 2 ) : k 1 , k 2 ∈ N ij ∧ (k 1 , k 2 ) ∈ E}.
That is, G ij is the subgraph that includes all the nodes (with their corresponding edges) whose indices are between i and j (both excluded). The number

of subgraphs G ij in G is (n -1) (n -2) 2
, the number of possible pairs (i, j). The maximum clique problem can be then solved as a maximization problem over the maximum cliques of all G ij . Let ÛG be the maximum clique in G. We can write:

Û = max ∀i,j:i<j ÛGij (1)
This decomposition does not seem, in principle, advantageous for efficiently computing the maximum clique in G. However, as we will see later, the computation of all ÛGij can be performed incrementally in such a way that the computation of a specific ÛGij can benefit from results previously calculated in terms of search space reduction.

To describe the incremental computation of the maximum cliques for a specific subgraph G ij let's assume that all values ÛG ik : i < k < j and ÛG kj : i < k < j have been previously computed. Let's also define the value B i,j as the number of Hamiltonian paths of length 3 (triangles) containing both i and j, and other intermediate nodes. We can then define an upper bound b i,j for ÛGi,j as:

b i,j = min( max ∀k:i<k<j∧(i,k)∈Eij ÛG i,k + 1, max ∀k:i<k<j∧(k,j)∈Eij ÛG k,j + 1, B i,j ) (2)
The meaning of Eq (2) can be explained as follows. The maximum clique size for G ij must be, at most, the size of the maximum clique for some G ik where i < k < j plus one, corresponding to the addition of the node j. Similarly, this size is limited by the size of the max clique in G kj and the addition of the node i. The most restrictive condition has to be satisfied and, therefore, the minimum of these two values is taken. Finally, there is a third condition given by the value B ij , which is also an upper bound for | ÛGij |. The inclusion of this third restriction will be discussed later as it constitute an initial point to compute the clique sizes for all the G ij .

From the previous discussion, we can also notice a necessary condition for a node k (such that i < k < j) to be a member of a clique ÛGi,j of size b, which is:

| ÛG i,k | + | ÛG k,j | ≥ b -1 (3)
This condition is trivial to prove because, if there is a clique ÛGi,j of size s, where k participates in that clique, then it is necessary that the sum of the sizes of the maximum clique from i to k and from k to j plus the addition of the node k itself (note that k is not included in either G i,k or Gk, j) is greater or equal to b. Hence, we can construct a subgraph of G i,j , which we call HV subgraph, composed of those nodes in G ij that fulfill the condition in Eq [START_REF] Bomze | The Maximum Clique Problem[END_REF], where b is approximated by the upper bound b i,j as it is defined in Eq [START_REF] Balas | Finding a maximum clique in an arbitrary graph[END_REF].

Thus, we can build the HV subgraph as follows:

HV i,j,b = {k : i < k < j ∧ (i, k) ∈ E ij ∧ (k, j) ∈ E ij ∧ | ÛG ik | + | ÛG k,j | ≥ b -1} (4)
Notice that the HV -subgraph for a specific G ij , given the bound b, consists of those nodes k in G ij that can form a clique of size b according to the clique sizes in previous G ik and G kj .From this, we can use this HV -subgraph to find the real value of | ÛGij | using any available algorithm for solving the maximum clique problem. The main advantage of using the HV graph here is that we can focus on those nodes that can actually contribute to a clique of size b thereby reducing the search space during the combinatorial search. In addition, we know that we have to finish the search as soon as we find a clique of size b. We can also notice that, by adopting this approach, apart from obtaining the maximum clique in G as it is stated in we get as a b product all the intermediate results ÛGij .

Example. Let G be a complete, 5-node graph where N = {1, 2, 3, 4, 5}. Then,| ÛG1,5 | = 3, because there are 3 intermediate nodes between 1 and 5.

Thus, HV 1,5,3 = {2, 3, 4}, because | ÛG 1,k | + | ÛG k,5 | ≥ 2
In the following sections, specific details about the algorithm for finding the maximum clique based on this framework along with some implementation issues are given.

An algorithm for finding the maximum clique

In this section we are going to rely on the framework presented in Section 3.1 to devise an algorithm for finding the maximum clique in a graph. We will start by giving some basic notation. Given a graph G = (N, E), we define the adjacency matrix A as a bidimensional matrix where each element (i, j) in A takes a value of 1 if (i, j) ∈ E and 0 otherwise. Without loss of generality and assuming a non-directed graph, we can consider A a triangular matrix. This way, we will assume that A i,j = 0 when j < i.

Next, we can define a new matrix, B obtained as B = A 2 &A where & represents the bitwise and operator that is an entrywise operator defined as:

(M 1&M 2) i,j = M 1 i,j If M 2 i,j = 0 0 otherwise (5)
And, therefore, the B matrix is computed as:

B i,j = A 2 i,j If A i,j = 0 ∧ i < j 0 otherwise (6)
The use of a triangular matrix is due to the fact that the element A 2 ij represents the number of paths of length 2 between the nodes i and j, and with the additional condition of existence of an edge (i, j), the element B ij represents the number of triangles containing that edge. The rationale behind the computation of this matrix B is that a clique of size n that has an edge (i, j), where i < j, also contains a sequence of n -2 triangles, n -3 squares, n -4 pentagons, and so on, all of them with paths in ascendant order, that is, ascendant Hamiltonian paths. Now, we will define a third matrix, C that will store the results for the sub-problems that are being solved during the search for the maximum clique. This way, each element C ij will take the value | ÛGij |. Before describing in detail the way in which the C matrix is constructed, we can enumerate two simple properties:

1. The C Matrix is a triangular matrix. Since each element (i, j) in C stores the size of the maximum clique for the subgraph G ij , all the elements C i,j where j ≥ i will be equal to zero.

2. The elements in the diagonal C i,j where i = j -2 will be either 0 if (i, j) / ∈ E or B i,j otherwise.

The C matrix is constructed incrementally. Initially, we can fill the positions corresponding to the diagonal (j -2, j) from their respective values in the B matrix. Next, we will compute the elements for all the columns j = 1 . . . n in C in ascending row order. That is, C j-3,j , C j-4,j , . . . , C 1,j . It is necessary to strictly follow this order, since to compute the current element C ij we have to first estimate an upper bound for | Û Gij |, which is computed as it is shown in Eq 2, which requires that all the values in the right part of the equation have been already calculated b i,j = min( max

∀k:i<k<j∧A i,k =1∧A k,j =1 C i,k + 1, max ∀k:i<k<j:A k,j =1∧A k,j =1 C k,j + 1, B i,j ) (7)
As we can see, Equation ( 7) is equivalent to Equation 2 but introducing the C matrix. Notice that the ÛG i,k and ÛG k,j in Equation (2) are taken from row i and column j in this C matrix. Next, we can build the HV -subgraph from the values stored in the C matrix:

HV i,j,bi,j = {k : i < k < j ∧ A i,k = 1 ∧ A k,j = 1 ∧ C i.k + C k,j ≥ b i,j -1} (8)
Again, the HV subgraph is obtained from the row i and column j in the C matrix (HV stands for Horizontal-Vertical). According to Equation ( 8), apart from the obvious fact that only the nodes connected to i and j are considered, we have two terms extracted from the already computed part of C. The first term (C i,k in the equation) corresponds to the H (Horizontal) contribution of the C matrix and represents the maximum clique between the i node and the considered node k. The second term (C k,j ) corresponds to the V (vertical) contribution of the C matrix and denotes the maximum clique between the considered node k and the j node. This way, only those nodes whose H and V contributions sum at least b i,j -1 are candidates to participate in a clique of size b i,j as it is shown in Figure 1. 0 0 0 0 1 1 1 2 3 2 1 2 3 4 5 6 7 8 9 10 1 0 0 1 2 3 3 4 4 5 4 2 3 4 0 0 0 0 0 0 1 1 2 1 5 0 0 0 0 0 0 1 1 2 0 6 0 0 0 0 0 0 0 0 1 1 7 0 0 0 0 0 0 0 0 0 0 8 0 0 0 0 0 0 0 0 0 0 9 0 0 0 0 0 0 0 0 0 0 10 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 3 3 2 3 Fig. 1: Computation of the upper bound and the HV -subgraph. In this figure, we can see how both the upper bound and the HV subgraph for the position (2, 10) are computed. First, we can obtain the bound as the minimum of the maximums of the row 2 and column 10 plus 1 (as indicated in Eq (2)). We can assume, in this example that B10,2 ≥ 3 to consider all the elements in that Equation. Once that the bound has been obtained, the HV -subgraph consists of all the nodes whose sum of their H and V components is equal or greater than b -1 (where b = 3 in this example). In the figure. the H and V component for each node are graphically connected. For instance, the H component of the node 7 is found in the position (2, 7) while the V component for this node is found in the position [START_REF] Erdős | On the evolution of random graphs[END_REF][START_REF] Konc | An improved branch and bound algorithm for the maximum clique problem[END_REF]. This way HV2,10,3 = {3, 4, 6, 7, 8, 9}. Notice that the node 5 is not included in the HV subgraph since H(5)+V (5) = 1 < 2

From the upper bound and the HV -subgraph we can compute the real value of C ij . Nonetheless, we have to remark that there is a circular dependence between the bound b i,j and the HV -subgraph. This is because the HV -subgraph depends on the previously computed bound and, similarly, the nodes in the HV graph condition the presence or not of a clique of size b i,j (for instance, if the number of nodes in the HV -subgraph is lower than b i,j , we can update the bound to a lower value and, correspondingly, the HV -subgraph has to be rebuilt according to the new bound).

On account of this, finding the exact value for C i,j ( ÛGi,j ) entails performing an iterative process until a clique of size b i,j is found in the current subgraph. Therefore, assuming that we have an algorithm called kClique that can determine whether or not there is a clique of k nodes in the HV -subgraph, we can follow the Algorithm 1. end if 9: end while 10: C i,j = b Notice that step 3 corresponds to the k-clique problem. Owing to the fact that this problem is similar in time complexity to the maximum clique problem [START_REF] Vassilevska | Efficient algorithms for clique problems[END_REF] we can obtain as a byproduct (one of the) maximum cliques between i and j in this step. Once the C matrix is completely filled, we can easily find the maximum clique in the graph by returning the clique corresponding to the position with the highest value in the C matrix. It is trivial to observe that it is not necessary to store the cliques for all the positions in C. Instead, we just have to save the maximum clique found so far and after finishing the computation of all the elements in C, this will be the maximum clique in the graph (as it is stated in Eq (3.1).

Regarding the search algorithm for finding the the maximum clique given the current HV -subgraph, we can follow a straightforward approach by relying on a state-of-the-art algorithm. Let's say that we have a generic algorithm to solve the maximum clique problem (we will call this algorithm MaxClique from now on). This leads to an interesting situation because the HV method can take advantage of any improvements in the state of the art in the clique problem just by using a new MaxClique.

Besides, we can easily construct an experimental framework to test the general performance of the HV algorithm. The point here is that we can run the MaxClique algorithm on the input graph as a baseline scenario and then the HV method using the same MaxClique to find the maximum clique in each HVgraph.

Pruning the search. An alternative approach

The technique described above computes all the ÛGij to solve the maximization in Eq (3.1) (we will call this method ExactHV). Although the computation of those intermediate values can be useful in some scenarios, to obtain the maximum clique of a given graph, we could improve the algorithm performance if we just compute the real values for those G ij whose initial bound b i,j is greater than the maximum clique found so far (we can call this new approach PrunedHV). In other words and if the current bound b i,j for an element C ij is greater than the best clique found, we will call to MaxClique to determine whether or not there is a clique of this size. Otherwise we just take the current bound as the value for C ij as it is shown in Algorithm 2 end if 8: end while 9: C i,j = b At this point, it can be necessary to remark that using the initial bounds in some positions in C has no impact on the final result of the search. This is because using an upper bound in any C position instead of the true value only causes that the HV -subgraphs for future elements in C contain more nodes than actually needed to compute the maximum clique for this position, but the real value of the maximum clique remains the same.

By introducing this pruned search, we can reduce the number of calls to MaxClique at the expense of increasing the size of the HV -subgraph in each call. Because of the exponential complexity of the available MaxClique algorithms, we could expect that below certain threshold (defined in terms of a combination of graph size and density) PrunedHV is faster than ExactHV. Conversely, beyond this threshold ExactHV would be more efficient since the reduction in the subproblems size will have a greater impact on the performance than the higher number of calls to MaxClique. Experiments with both approaches (pruned and not pruned) will be reported in Section4.

In Algorithm 3 the full algorithm for finding the maximum clique in both versions (pruned and exact) is defined.

Using a cache of cliques

In the exactHV method, to compute each element in the C matrix we obtain as a byproduct a clique of the corresponding subgraph G ij . We can benefit from this by storing these found cliques in order to reduce the number of calls to MaxClique for future elements in C. This is based on the assumption that when computing the element C ij we can obtain a clique of size b = b ij by adding the nodes i and j to some of the cliques found previously. Specifically, given the current subgraph HV i, j, we can take those elements k ∈ HV i, j : H(k) ≥ b -1. That value b -1 of the H component indicates that we previously found a clique of size b -1 in the element C i,k . If this clique was stored, we can try to get a bigger clique (of size b) by adding the node j to this previously saved clique. We can follow a similar reasoning for those elements l ∈ HV i, j : V (l) ≥ b -1 in such a way that we know that a clique of size b -1 was found in the element C l,j and, therefore, we can add the node i to that clique to try to obtain a clique of size b.

Notice that to estimate the correct value of each C ij element it is enough to find a clique of the size indicated by the current bound. As a result of this, it seems plausible that we can save some calls to MaxClique by trying to extend those cliques found previously. In Algorithm 4 the use of the cache of cliques is explained.

Algorithm 3 Algorithm for finding the maximum clique. The algorithm takes as inputs the graph G, a generic method MaxClique and a parameter that indicates which one of the two algorithms (ExactHV or PrunedHV) will be carried out HV = HV i,j,b as stated in Eq (4)

13:

C i,j = b 14:
if prune=true then Update HV = HV i,j,b as stated in Eq (4) return bestClique 43: end function Algorithm 4 Using a cache of cliques to reduce the number of calls to MaxClique. This algorithm would be inserted just before the calls to Maxclique in Algorithm 3 in lines 16 and 25.

1: for each node k ∈ HV i,j,bi,j :

H[k] ≥ b(i, j) -1 do 2: candidate → clique stored in position (i,k) 3: candidate → candidate ∪ j 4:
if candidate is a clique then 5:

C(i, j) = b(i, j) return 6:
end if 7: end for 8: for each node l ∈ HV i,j,bi,j :

V [l] ≥ b(i, j) -1 do 9:
candidate=clique stored in position (l,j)

C(i, j) = b(i, j) return 13:
end if 14: end for

Experimental Framework

In this section, we will describe the experimental framework adopted in order to test the efficiency of the approach presented here. For these experiments, we have used Tomita's MCS [START_REF] Tomita | A simple and faster branch-and-bound algorithm for finding a maximum clique[END_REF] as the baseline scenario. Besides, the MaxClique algorithm used to compute the real value of each ÛGij (C ij ) as it is shown in lines 16 and 25 in Algorithm 3 is, as well, MCS. The experiments have been performed on two different sets. First, several random graphs with different sizes and densities were produced. Second, a subset of the DIMACS benchmark set was employed as well.

In Table 1 the main results on random graphs, generated according to the well known Erdős-Rényi model [START_REF] Erdős | On random graphs i[END_REF], [START_REF] Erdős | On the evolution of random graphs[END_REF], are summarized. Each entry in the table represents the average time for 10 instances of each graph.

From the results presented in Table 1, we can conclude that the HV approach is not appropriate for small graphs and low densities (specifically, M CS is always faster for graphs with 200 or fewer nodes). This can be explained by the fact that the polynomial cost associated to computing the the C matrix elements and obtaining the HV -subgraphs, which is O(n 3 ) heavily influences the performance of the algorithm in this cases. In addition, the decomposition in subproblems can be not so advantageous for this small size problems. However, as the problem complexity increases (either in size or density) we can see how the difference in performance between M CS and HV becomes tighter and tighter to reach a point where HV outperforms M CS.

In spite of the fact that the improvements presented in Table 1 can be considered modest, it is worth mentioning that as the size or density of the input random graphs grow, we can see how the performance gap between M CS and HV also grows as it is shown in Figure 2. In this figure, we can see how the the time complexity for HV grows at a lower rate than M CS when either the size or graph density is increased. Therefore, we can expect notable improvements in performance for very large graphs (or very dense graphs given a minimum number of nodes). In this case, the decomposition in smaller sub-problems is clearly beneficial given the exponential nature of the available algorithms for the clique problem. In Table 2 we can see the number of size of subproblems for several random graphs.

Tab. 1: Results on random graphs. PrunedHV corresponds to the technique described in Section 3.3 Regarding the results on the DIMACS benchmark, reported in Table 3 they show that, for some families of graphs, we can see a behavior similar to random graphs. Conversely, for some other graph types the HV algorithms shows poor performance. Families P-hat and Dsjc are examples of this and,in particular, Mann graphs. In general, we can say that for those graphs that have stronger symmetries, the difference is notable in favor of MCS. The only exception is the Johnson16-2-4 graph which, despite being a small graph, is surprisingly solved faster by HV.

Graph

MCS

From this results, if we consider the Erdős-Rényi model as a possible distribution of all the different classes of graphs, the HV algorithm is a new interesting approach to cope with the problem of the maximum clique, even considering the poor performance showed for some graph types. Moreover, this approach can be considered as a kind of meta-algorithm since any method to solve the clique problem can be embedded into any of the HV variants.

Finally, we have to mention that, from the algorithm presented in Algorithm 3 an interesting question arises. The calls to MaxClique (which are calls to M CS in the experiments here reported) can be replaced with recursive calls to the HV algorithm itself. Although we can not expect any improvement for the graph sizes and densities considered in this section, we can guess that, for very large graphs, the inclusion of these recursive calls can improve the HV results shown in Table 1. This improvement can be even higher if the base case for recursion is defined as a call to other MaxClique algorithm (for instance, M CS) when the size and density of the current HV -subgraph is below a fix threshold.

Complexity

As commented in the introduction section, the algorithmic complexity for the worst case for maximum clique searching algorithms is exponential, but no explanation is due, in general, about the average complexity because the difficulty to cope with this case. Thus, no analytical studies showing the complexity of this problem can be found in the literature. On the contrary, our approach allows us to achieve an approximation, albeit still preliminary and not very developed, about this aspect.

To do this, we consider as hypothesis a distribution of cliques on graphs following the well known Erdős-Rényi model [START_REF] Erdős | On random graphs i[END_REF], [START_REF] Erdős | On the evolution of random graphs[END_REF]. According to this model, every edge between two nodes will be present with a certain probability p. We consider the hypothesis of p = 0.5 to calculate the time complexity under this assumption.

We call T HV (n) the time required to compute the solution with HV algorithm for a problem of size n in that hypothesis, where this value is the number of nodes in the graph. Of course, if this number is small, the complexity will be constant, in O(1). Otherwise, the first part of the algorithm consists on a matrix multiplication, and this operation can be carried out in a time O(n 3 ) or smaller, thus, in the worst case we use this value as an upper bound of the complexity. After that, we take the calculations of the values for each C(i, j). As there are a number n • (n -1) of different pairs (i, j), that is the number of calls needed for the algorithm to compute those values C(i, j) using the MCS algorithm. We call T M CS (n) its complexity, but now, the size of the problem with which we feed this algorithm becomes n/K for a certain value K. This value is unknown, but we can use the hypothesis of the Erdős-Rényi model with p = 0.5, and, in this case, we can guarantee that this number is greater than 4, because this is the expected value of number of nodes with connection with both nodes i and j.

Thus, we have this equation to approximate the time complexity:

T HV (n) ∈ O(n 3 ) + n • (n -1) • T M CS (n/K) If n > Size O(1)

Otherwise

The value of T M CS (n) is also unknown, but as we rely on the MCS algorithm to solve the maximum clique problem between nodes i, j, in the same way, we can consider a recursive version of HV , using the HV approach itself to solve the corresponding subproblems. We can calculate the complexity of this new algorithm by means of a differences equation:

T HV (n) ∈ O(n 3 ) + n • (n -1) • T HV (n/K) If n > Size O(1)

Otherwise

To solve this recurrence equation, we assume, in accordance with the considerations indicated above, that K > 4. In fact, the value, depending of n, is greater than 4, but, again, we can consider this value as upper bound, and the solution of this equation in that case is

T HV (n) ∈ n O(log n)
Finally, experimentally, we can check that time performance of Tomita algorithm is better than HV recursive algorithm, and thus we can consider

T T (n) < T HV (n) ∈ n O(log n)
Of course, this calculation is with respect to Erdős-Rényi graphs, and the hypothesis p = 0.5, although the result is valid for any other probability, taking into the account that this variation will influence only in the base of the logarithm.

Conclusions and Future Work

In this work, we have presented a new technique to solve the classical problem of finding the maximum clique in a graph, which is a well-known N P -complete problem. This algorithm is quite novel, since the traditional successful methods are based on backtracking and pruning methods, whereas the method here proposed is based on a heterodox version of the divide and conquer methodology.

As a result of decomposing the main problem into smaller elements, we can obtain, in the first place, additional information in the form of sub-cliques for certain subgraphs of the main graph. In addition and, more interestingly, the algorithm seems to outperforms the classical backtracking approaches for random graphs where the size or density are high enough. Although the improvements obtained are not actually notable we can expect that, asymptotically, the difference in performance is quite remarkable.

Regarding future work, experiments with larger graphs have to be conducted in order to reinforce the conclusions presented about the algorithm performance. In addition, the results for some DIMACS graph have to be analyzed to find possible improvements when dealing with graphs having strong regularities. Finally, a deep theoretical and empirical study of the recursive version of the HV algorithm, as described in Section 4 could be addressed as well.

Algorithm 1

 1 Computing the final value of C ij 1: Compute an initial bound b = b i,j as it is stated in Equation (7

Algorithm 2

 2 Computing the final value of C ij . Pruned version 1: Compute an initial bound b = b i,j as it is stated in Equation (7). 2: res=0 3: while b < bestCliqueF ound and res < b do

1 :while i ≥ 1 do 10 :

 110 function MAXIMUMHVCLIQUE(G,maxClique,prune) if B i,j > 0 then 11: b = b ij as stated in Eq (2) 12:

  HV = HV i,j,b as stated in Eq (

  Size and number of subproblems for several random graphs. For PrunedHV (HVMax) and ExactHV(HVMax) The number of subproblems (calls to MaxClique) along with the maximum and average subproblem size (the subproblem size corresponds to the size of the subgraphs in the calls to MaxClique) are shown In this figure we can see a comparison between the performance of MCS and HVApprox for random graphs as the graph size increases. The left figure shows the MCS and HVApprox times for different graph sizes generated with p = 0.5 while the right figure shows the performance improvement in percentage. Results on the DIMACS benchmark. HVMax corresponds to the technique described in Section 3.3. The blank entries indicate that the algorithm was stopped because the time elapsed was too high

	PrunedHV ExactHV 2.8 7.0 2.8 7.0 30.6 99.9 11.8 54.4 0.5 15 -16 234.2 182.5 p w 0.4 10-12 2.2 0.5 13 2.2 0.6 16-17 33.2 1000 0.4 12 n 500 10.6 591.6 0.6 19-20 12295 10635 28401 1500 0.4 12-13 127.0 118.2 279.4 0.5 16 -17 4930 3674 10689 2000 0.4 13 860.9 724.9 1575 0.5 17 50151 34188 101806 3000 0.4 14 12512 9233 21020 PrunedHV (HVMax) ExactHV (HVMax) #subpr. max size av. size #subpr. max size av. size 0.85 204190 p 177 107 228954 177 96 0.7 975751 184 90.5 1062611 184 84 1000 0.6 Tab. 2: n 300 500 7841222 228 103 10005828 228 94 0 10000 20000 30000 40000 50000 400 600 800 1000 1200 1400 1600 1800 2000 time (s) n MCS HVApprox -30 -20 -10 0 10 20 30 40 400 600 800 1000 1200 1400 1600 1800 Improvement (%) n HVApprox/MCS MCS PrunedHV ExactHV Instance n w brock400 1 400 27 604.7 399.3 2384 brock400 2 400 29 272.2 352.9 1432 brock400 3 400 31 415.1 157.8 1328 brock800 1 800 23 7916 5937 18812 brock800 2 800 24 7216 8873 21474 brock800 3 800 25 4860 3671 18588 brock800 4 800 26 3669 5041 19279 phat500-3 500 50 137.9 319.5 378.38 phat700-2 700 44 4.4 15.07 38.3 phat1000-2 1000 46 197.7 440 1034 phat1500-1 1500 12 3.1 10.3 27.9 Fig. 2: Tab. 3: Graph C250.9 250 44 3827 3761 9146	2000

| ÛG1,2 | = 0 | ÛG2,5 | = 2 | ÛG1,3 | = 1 | ÛG3,5 | = 1 | ÛG1,4 | = 2 | ÛG4,5 | = 0

Acknowledgements

This work has been partially supported by Spanish Minister of Economy project TIN2015-65845-C3-2-R.


Shuda Yu

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

Incremental Reconstruction of Manifold Surface from Sparse Visual Mapping

Automatic image-based-modeling usually has two steps: Structure from Motion (SfM) and the estimation of a triangulated surface. The former provides camera poses and a sparse point cloud. The latter usually involves dense stereo. From the computational standpoint, it would be nice to avoid dense stereo and estimate the surface from the sparse cloud directly. Furthermore, it would be useful for online applications to update the surface while the camera is moving in the scene. This paper deals with both requirements: it introduces an incremental method which reconstructs a surface from a sparse cloud estimated by incremental SfM. The context is new and difficult since we ensure the resulting surface to be manifold at all times. The manifold property is important since it is needed by differential operators involved in surface refinements. We have experimented with a hand-held omnidirectional camera moving in a city.

Introduction

The estimation of the scene surface viewed by a camera is an important requirement for applications such as augmented reality, collision detection while moving along a planned path, environment modeling etc. However, the majority of methods which produce a 2-manifold (e.g. dense multi-view stereo [START_REF] Seitz | A comparison and evaluation of multi-view stereo reconstruction algorithms[END_REF] or methods based on [START_REF] Kazhdan | Poisson surface reconstruction[END_REF]) are batch methods. They are not easy to adapt to the incremental context where a surface is computed at all times from the progressive availability of images or 3d points.

In a 2-manifold (i.e. 2d topological manifold surface), each point of the surface has a surface neighborhood which is homeomorphic to a disk. Thus, a triangulated 2-manifold is not a simple soup of badly connected triangles. Each triangle is exactly connected by its 3 edges to 3 other triangles, the surface has neither holes nor self-intersections, and it cuts R 3 into free-space and matter regions. Here we assume that the 2-manifold is triangulated.

The manifold property provides the possibility of estimating safely the differential operators [START_REF] Meyer | Discrete differential-geometry operators for triangulated 2-manifolds[END_REF] used by surface refinements like surface fairing and dense stereo. In this paper, we use a surface denoising based on discrete Laplacian [START_REF] Taubin | A signal processing approach to fair surface design[END_REF], whose performances are degraded if we apply it on a non-manifold surface. Fig. 1 shows an example for car and street (during winter) reconstructed by our methods from a sparse SfM point cloud. More generally, a lot of Computer Graphic algorithms do not apply if the triangle list is not a 2-manifold [START_REF] Botsch | Polygon Mesh Processing[END_REF]. Note that a 2-manifold directly estimated from the sparse cloud produced by SfM would be ideal for both time and space complexities. This surface would also be useful for initializing a more accurate (but more costly) surface reconstruction method such as surface deformation minimizing photo-consistency [START_REF] Hiep | Towards high-resolution large-scale multi-view stereo[END_REF][START_REF] Newcombe | Live Dense Reconstruction with a Single Moving Camera[END_REF]. The dense-stereo/deformation step is outside the scope of this paper.

To our knowledge, this paper presents the first incremental method which provides a 2-manifold from a sparse cloud of reconstructed interest points provided by SfM. Here "incremental" means that a 2-manifold obtained before time t is locally updated using 3d points provided by SfM at t to obtain the 2-manifold at t. SfM is also incremental: the geometry at t (camera poses and sparse point cloud of the sequence up to t) is a local update of a geometry before t. We focus on the mapping scenario such as a camera mounted on vehicle/robot/human exploring an unknown and large scene. In our experiments, a hand-held omnidirec- tional (calibrated) video camera is moving in a city. This is different to the scenario of a camera moving in a limited workspace such as desk-like/indoor environments [START_REF] Newcombe | Live Dense Reconstruction with a Single Moving Camera[END_REF].

1 (A) (B) (C) F F F F F F F F F F F F F F F F F F F F F F (D)
Section 2 compares our work against others. Sections 3 and 4 describe a batch method [START_REF] Yu | Surface reconstruction of scenes using a catadiopric camera[END_REF] and our incremental method, respectively. The former is useful to describe the latter. Lastly, the experiments and conclusion are given in Sections 5 and 6.

Previous Work

We discuss the previous works which reconstruct a surface from sparse (not dense) cloud of features reconstructed by SfM. Since most of these methods use a 3d Delaunay triangulation, we start with a reminder of what this is.

A lot of surface reconstruction methods [START_REF] Cazals | Delaunay triangulation based surface reconstruction: ideas and algorithms[END_REF] are based on the following property [START_REF] Amenta | Surface reconstruction by voronoi filtering[END_REF]: if P is a sufficiently dense sample of points on the (smooth) scene surface, a good approximation of this surface is given by a list of triangles of the 3d Delaunay triangulation T of P . T is a list of tetrahedra which partitions the convex hull of P such that (1) the vertices of all tetrahedra are P and (2) the tetrahedra circumspheres do not contain a vertex within them. The triangles are the facets of the tetrahedra. Fig. 2 (A andB) illustrates this property in the 2d case.

In our case, the points are reconstructed from images and we know the list V i of indices of images which reconstruct the 3d point p i ∈ P . We refer to a ray as a line segment linking p i to the j-th camera location c j such that j ∈ V i .

The free-space carving methods [START_REF] Faugeras | Representing stereo data with the delaunay triangulation[END_REF][START_REF] Labatut | Efficient multi-view reconstruction of large-scale scenes using interest points, Delaunay triangulation and Graph-Cuts[END_REF][START_REF] Pan | ProFORMA: probabilistic feature-based on-line rapid model acquisition[END_REF][START_REF] Lovi | Incremental free-space carving for real-time 3d reconstruction[END_REF][START_REF] Yu | Surface reconstruction of scenes using a catadiopric camera[END_REF] use rays as visibility constraints to label the tetrahedra of a 3d Delaunay: a tetrahedron is free-space if it is intersected by at least one ray, otherwise it is matter. These methods deal with sparse cloud of points (or edges [START_REF] Faugeras | Representing stereo data with the delaunay triangulation[END_REF]), and only [START_REF] Lovi | Incremental free-space carving for real-time 3d reconstruction[END_REF] is incremental. They (except [START_REF] Faugeras | Representing stereo data with the delaunay triangulation[END_REF] and [START_REF] Yu | Surface reconstruction of scenes using a catadiopric camera[END_REF]) directly consider the surface as the list of triangles separating the free-space and matter tetrahedra (see Fig. 2). Unfortunately, the resulting surface may be non-manifold. For example, the surface has a singularity at vertex v if all tetrahedra which have vertex v are matter, except two free-space tetrahedra ∆ 1 and ∆ 2 such that the intersection of ∆ 1 and ∆ 2 is exactly v. In [START_REF] Faugeras | Representing stereo data with the delaunay triangulation[END_REF] ( [START_REF] Yu | Surface reconstruction of scenes using a catadiopric camera[END_REF], respectively), a region growing procedure in the list of matter (free-space, respectively) tetrahedra removes all singularities and provides a 2-manifold.

Methods [START_REF] Morris | Image-consistent surface triangulation[END_REF][START_REF] Taylor | Surface reconstruction from feature based stereo[END_REF][START_REF] Hilton | Scene modelling from sparse 3d data[END_REF] use 2d Delaunay triangulations and deal with a sparse cloud of features. Only [START_REF] Hilton | Scene modelling from sparse 3d data[END_REF] is incremental, but the surface is not manifold (it may have holes) and the approach is applied to a small sequence of real images.

Batch Surface Reconstruction

The batch method is useful to describe our incremental method. Its steps are 3d Delaunay, Ray tracing, 2-Manifold Extraction, Topology Extension, and Surface Denoising.

3d Delaunay Triangulation T Assume that SfM estimates the geometry of the whole image sequence. The geometry includes the sparse cloud of points {p i }, camera locations {c j } and rays defined by visibility lists {V i } (notations in Section 2). Point p i has poor accuracy if it is reconstructed in degenerate configuration [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]: if p i and all c j , j ∈ V i are nearly collinear. This case occurs in part of the camera trajectory which is a straight line and if points reconstructed from this part are close to the straight line. Thus, p i is added in T if and only if there is an angle c j p i c k (j, k ∈ V i ) larger than threshold ǫ.

We also add extra points in T . For the clarity of our current paper, we only mention that (1) these extra points act reconstructed points with empty visibility lists and (2) they are randomly added in the neighborhood of the camera trajectory. The reason (removal of spurious arks/handles) and technical details are in Section 4 of our previous paper [START_REF] Yu | Surface reconstruction of scenes using a catadiopric camera[END_REF].

Ray Tracing As all tetrahedra are initialized matter, raytracing is applied to each ray to force into free-space all tetrahedra intersected by the ray. T is defined by a graph: a graph vertex is a tetrahedron, a graph edge is a triangle between two tetrahedra. Tracing a ray c j p i is a walk in the graph, starting from the tetrahedron which contains c j , moving to another tetrahedron through the triangle intersected by the line segment c j p i , stopping in tetrahedron which has vertex p i . Now we know the label (matter or free-space) of all tetrahedra which partition the convex hull C of the points, but the label of R 3 \ C is unknown. In our case, the points are reconstructed in almost all directions around view points (we reconstruct an environment). Thus, Figure 3. v is regular since the edges opposite to v define a simple polygon abcdef ga on the surface. v ′ and v ′′ are not regular since polygons a

′ b ′ c ′ d ′ e ′ f ′ g ′ a ′ -p ′ q ′ r ′ t ′ p ′ and a ′′ b ′′ c ′′ d ′′ e ′′ f ′′ g ′′ p ′′ q ′′ g
′′ a ′′ are not simple (the former is not connected, the latter has multiple vertex g ′′ ).

view points and rays are in the convex hull of the points. Since the rays do not intersect R 3 \ C, R 3 \ C is matter.

2-Manifold Test

The target surface S is a list of triangles of T which should be 2-manifold. Let v be a point in S. We say that v is regular if it has a neighborhood in S which is topologically a disk. Otherwise v is singular. By definition, S is 2-manifold if all its points are regular. In our context where S is a list of triangles of T , it is sufficient to check that each vertex v of S is regular using the following neighborhood of v: the list of the S triangles which have vertex v [START_REF] Boissonnat | Geometric structures for three-dimensional shape representation[END_REF]. Then, v is regular if and only if the edges opposite to v in the triangles of S having v as vertex form a simple polygon (Fig. 3). A simple polygon is topologically a circle, i.e. a list of segments which forms a closed path without self-intersection.

2-Manifold Extraction

The target surface S should also separate free-space and matter as far as possible, under the constraint that it is 2-manifold. A 2-manifold cuts R 3 in regions labeled outside (outside the matter) and inside. Here the outside region O contains a maximum of free-space tetrahedra and does not contain matter tetrahedron. A region growing process is used: O grows from ∅ by adding free-space tetrahedra one by one, such that the border δO of O remains 2-manifold. The final δO is S. We know that δO is 2-manifold if and only if all δO vertices are regular. Since tetrahedra are added one by one, the neighborhoods of at most four vertices of δO (those of the added tetrahedron) are modified. So we only need to check that these vertices are regular after the tetrahedron is added. If this is not the case, this tetrahedron is removed from O and we try another one. The final O depends on the addition order of the tetrahedra in O. We choose the added free-space tetrahedron such that it has a facet included in δO, i.e. it is in the neighborhood of O. A priority is also defined for each free-space tetrahedron: the number of rays which intersect the tetrahedron. The tetrahedra in the neighborhood of O are stored in a heap (priority queue) for fast selection of the tetrahedron with largest priority.

The inputs of the region growing are the initial O, the set F 0 of tetrahedra where the growing is possible, the set Q 0 of tetrahedra which includes the initial value of the heap, and function r which maps a tetrahedron ∆ to the number of rays which intersect ∆. In the batch case, O = ∅, Q 0 = T and F 0 is the list of free-space tetrahedra of T . The output is O. Here is the algorithm in C style.

if (∆ ′ ∈ F 0 and ∆ ′ / ∈ O) Q ← Q ∪ {∆ ′ }; } else O ← O \ {∆}; }
Topology Extension The δO genus can not be changed if the tetrahedra are added one by one by the algorithm above: O always has the ball topology. This is problematic if the true outside does not have the ball topology, e.g. if the camera trajectory contains closed loop(s) around building(s). In the simplest case of one loop, the true outside has the toroid topology and the computed outside O can not close the loop. This problem is corrected as follows. First, we find a vertex in δO such that all inside tetrahedra incident to this vertex are free-space. Second, we force all these tetrahedra to outside (O is increased). Third, we check that all vertices of these tetrahedra are regular. In case of failure(s), these tetrahedra are restored to inside (O is decreased). Finally, we alternate this scheme and the previous algorithm until no more tetrahedron can be added in O. Here, Q 0 is the list of tetrahedra neighbors of the forced tetrahedra above, and F 0 is unchanged.

Surface Denoising

The S reconstruction noise is reduced thanks to a smoothing filter p ′ = p + ∆p where p is a vertex of S and ∆p is a discrete Laplacian defined on S vertices [START_REF] Taubin | A signal processing approach to fair surface design[END_REF]. The smoothed p ′ is stored in a distinct array of p. We don't apply p ← p ′ to avoid the computation overhead due to vertex update in T .

Incremental Surface Reconstruction

Our method is defined by a main loop which alternates the incremental versions of the steps in Section 3. Integer t specifies the current time and the keyframe index.

Incremental SfM First, a new keyframe is selected from the input video and interest points are matched with the previous keyframe using correlation. The new keyframe is such that the number of its matches with the two previous keyframes is larger than a threshold. Then, the new pose is robustly estimated (using Grunert's method and RANSAC) and new 3d points are reconstructed from the new matches. Lastly, local bundle adjustment refines the geometry of the l-most recent keyframes. Using l = 3, the l most recent keyframes are t -2, t -1 and t. No more details are given on this SfM step since it is similar to [START_REF] Mouragnon | Generic and real-time structure from motion using local bundle adjustment[END_REF].

3d Delaunay Triangulation T We add point p to T once p reaches its final value by SfM. From the computational standpoint, this is efficient since we don't need to update T every time SfM updates p. More precisely, we add in T at time t (after the SfM step at t) every 3d point p such that its last track is in keyframe t -2. As in the batch case, we check that p is not in a degenerate configuration and we add extra points. A fixed number of extra points is added in the neighborhood of c t .

Dating Our incremental method needs a creation date for all tetrahedra and vertices. Furthermore, we use the Delaunay implementation of CGAL [1] which adds the points one by one to T . The adding of p destroys list L d (p) of tetrahedra and creates list L c (p) of other tetrahedra. So we assign date t to p and to the tetrahedra of L c (p) if p is added to T at time t. We also need the smallest date d t of all outside tetrahedra in L d (p) destroyed at time t (for all p added at t). Each tetrahedron was labeled outside or inside by the 2manifold step below. Both lists L c (p) and L d (p) are easy to compute thanks to CGAL functions.

Ray Tracing Tracing all rays available at date t is too time consuming. Here we do the following approximation: the label (free-space or matter) of a tetrahedron is defined by the rays which have creation dates similar to or greater than that of the tetrahedron. A ray has the creation date of its 3d point, defined in the "Dating" step. According to this approximation, we only need to ray-trace the most recent rays. At date t, we apply ray-tracing to the small list of rays which have creation dates in {tk, • • • , t -1, t}, where k is a threshold.

2-Manifold Extraction

Starting the region growing from O = ∅ as in the batch case is too time consuming. In the incremental case, we propose a method which starts the region growing from a list O obtained at a recent date. We regroup the free-space tetrahedra into different layers L t ′ by creation date t ′ and the idea of growing the outside O by layer of creation date comes naturally to mind. We grow O layer by layer, and for each layer, only free-space tetrahedra created before or at this layer can be added into O. As a result, for each layer, we could extract a 2-manifold as the border of the tetrahedra list O. Then the 2-manifold of the next layer can be easily computed by starting from that of the current layer. In practice, we prefer to grow by pack of n layers for efficiency. At each time t, the method holds several lists of outside tetrahedra which correspond to particular creation dates (multiples of n): O n , O 2n , ... O itn where i t is the largest integer such that i t n ≤ t. These lists and another list O t of outside tetrahedra meet

O n ⊆ O 2n ⊆ • • • ⊆ O (it-1)n ⊆ O itn ⊆ O t ∀t ′ ∈ {n, 2n, • • • , i t n, t}, O t ′ ⊆ L 1 ∪ L 2 ∪ • • • ∪ L t ′ and the border of O t ′ is 2-manifold.
(1)

The border of O t is the target 2-manifold S at time t. Furthermore, the region growing at time t + 1 is started from one of the O t ′ lists above.

To simplify notations in this Section, time t starts from 1 (not 0) and we define O 0 = ∅. If t ≤ n, we apply the batch region growing in all free-space tetrahedra from O 0 to obtain O t . Now assume that t > n. The algorithm works from Eq. 1 at t -1 to Eq. 1 at t. Fig. 4 illustrates this if t = 98, n = 20, d t = 51. Remember that the point additions at time t destroy tetrahedra and we know the smallest date d t of the destroyed outside tetrahedra. Let i 0 be the largest integer such that i 0 n < d t . If i ≤ i 0 , O in is unchanged and its border is still manifold. If i 0 < i, tetrahedra may be destroyed in O in , its border may be non manifold and O in should be recomputed. Time starts from 1, thus 1 ≤ d t , 0 ≤ i 0 and O i0n exists. Then we apply region growing (2-manifold extraction in Section 3) from O i0n to obtain O (i0+1)n . We also apply region growing from O (i0+1)n to obtain O (i0+2)n and so on, until we obtain O itn . Lastly, we apply region growing from O itn to obtain O t . Remember that F 0 and Q 0 should also be defined for region growing from O in to obtain O (i+1)n (or O t ), as mentioned in Sec-tion 3. For time complexity reason, we don't use

Q 0 = T but the most recent layers Q 0 = L in-b0 ∪ • • • ∪ L (i+1)n
where b 0 ∈ N is constant. We also define F 0 by the freespace tetrahedra of

L (i0-1)n ∪ • • • ∪ L (i+1)n (not those of L 0 ∪ • • • ∪ L (i+1)n ).

Topology Extension

The O in (including O t ) obtained at the previous step are improved by an incremental version of "Topology Extension" (Section 3) after the region growing from O (i-1)n to O in . The improved O in still meet Eq. 1. "Topology Extension" is only applied to the most recent vertices of S which have creation dates inb 1 , • • • , in -1, in where b 1 ∈ N is a threshold. These vertices are only tried once (there is only one "2-Manifold Extraction" and one "Topology Extension" for a given i). Here we use Q 0 of Section 3 and F 0 is the list of free-space tetrahedra of

L (i0-1)n ∪ • • • ∪ L in .
Surface Denoising Denoising all vertices of S is too time consuming. In the incremental case, we only need to smooth vertex p of S if its smoothing p ′ at time t is different to that at t -1 due to the steps above. As in Section 3, the smoothing is p ′ = p + ∆p where ∆p only depends on p and N (p). Neighborhood N (p) is the list of vertices which are connected to p by an edge of S. Thus p ′ is (re)calculated if p is a new vertex of S or if N (p) changes at t. The tetrahedra list O t \O i0n is the outside volume grown by steps "2-Manifold Extraction" and "Topology Extension" at t. Furthermore, all S changes at t are on the border of O t \ O i0n . So we (re)-calculate p ′ if N (p) ∪ {p} contains at least one vertex of the border of O t \ O i0n .

Experiments

Synthetic Sequence

Here we compare the performances of the batch (Section 3) and the incremental (Section 4) surface reconstruction methods on the same sparse cloud of 3d points estimated from images of a synthetic scene. The synthetic scene is manually generated from real images taken in a city. The trajectory is a 230 m long closed loop around a building including several shops. The images are generated by ray-tracing and taking into account the ray reflection on the mirror. The catadioptric camera has axial symmetry. The large circle, which contains the scene projection in the image, has a 600 pixel radius. Fig. 5 (top-left corner) shows two images of the synthetic sequence.

SfM [START_REF] Lhuillier | Automatic scene structure and camera motion using a catadioptric camera[END_REF] reconstructs 600 camera poses and a sparse cloud of 257336 3d points from the sequence. We approximate the true calibration by a central model and refine the radial distortion parameters using bundle adjustment. We estimate the similarity transformation R which minimizes

E(R) = 599 i=0 ||R(c i ) -c g i || 2
, where c i and c g i are the es- Qualitative Comparison Before using the ground truth, it is interesting to compare batch and incremental results. The 3d Delaunay triangulation has 123196 vertices and 750219 tetrahedra. The numbers of free-space tetrahedra are 494562 for batch and 490174 for incremental (the difference is 0.9%). The batch and incremental (final) surfaces have 233273 and 231826 triangles, respectively.

Remember that the list O of outside tetrahedra grows in the list of free-space tetrahedra. Thus the ratio between the numbers of outside and free-space tetrahedra can be used to compare the performances of the growing steps (2-Manifold Extraction and Topology Extension). Batch surface has 89.1% and incremental (final) surface has 85.8%. We explain this results as follows: the incremental growing is more constrained than the batch growing. In the incremental case, both dating and manifold constraints are used. The batch method only uses manifold constraints. In prac-Table 1. Errors of batch and incremental (final) surfaces using ground truth surface. The numbers between parentheses are obtained for twice smaller images. method inliers (%) median (cm) 90% quantile (cm) batch 75.7 (73.9) 8.0 (64.3) 55 (103) increm. 72.4 (69.7) 8.6 (71.0) 50 ( 106) tice, the ratio can not reach 100% since ray-tracing alone does not enforce the manifold constraint between free-space and matter tetrahedra.

Quantitative Comparison Now we define an error function to compare the estimated surface (batch or incremental) against the ground truth surface. At first glance, we could use distance e(p) between the ground truth surface and vertex p of the estimated surface [START_REF] Seitz | A comparison and evaluation of multi-view stereo reconstruction algorithms[END_REF]. Unfortunately, this error is biased in favor of reconstructed areas which have the largest densities of reconstructed points (ground parts have low textures and densities, walls have high densities). A second idea is the use of the same error such that p samples uniformly the estimated surface. However, this method has drawback since the closest point in the ground truth surface does not necessarily correspond to the same point p.

Our solution does not have the problems above. Let q be a pixel in an image of the sequence. Let p e be the intersection of the estimated surface and the back-projected ray of q by the estimated camera pose. Let p g be the intersection of the ground truth surface and the back-projected ray of q by the ground truth camera pose. In both cases, if there are several intersections, we take the intersection which is the closest to the camera pose. Then we use e(q) = ||p e -p g ||. If p g does not exist or e(q) > µ 0 (where µ 0 = 2 m), we assume that the point matching (p e , p g ) is outlier (e.g. for the pixels of the sky) and we ignore the error for q. In practice, we estimate the statistic of e(q) by uniform sampling of q in all images of the sequence. We sample 6000000 pixels in the sequence. Tab. 1 provides the results for both batch and incremental methods. We see that the batch method has slightly better results than the incremental method.

Lastly, the same experiment (both SfM and surface calculations) is re-done for the same images down-sampled by 2. We found E(R)/600 = 56 cm, which implies that the SfM drift is larger than in the previous case. According to Tab. 1, the batch surface is still the best and the surface accuracies are degraded. Fig. 5 shows the sparse point cloud by SfM and the surfaces.

Real Sequence

Our (equiangular) calibrated catadioptric camera is the 0-360 mirror mounted on the Canon Legria HSF10. We take a 1920*1080 AVCHD (MP4) video walking in a city during 505 seconds and pointing the mirror toward the sky by hand. Ground truth is not available, but we know that the trajectory length is about 800 m. The view field is 360 degrees in the horizontal plane and 51-58 degrees above and below. Fig. 6 shows our camera and several images of the sequence. The horizontal and vertical radii of the large ellipse, which contains the scene projection in the images, are 700 and 693 pixels, respectively.

The method in Section 4 is applied to the (down-sampled by 2) images with the same parameters as in Section 5.1. 1033 keyframes are selected from 25278 images. About 600 Harris points are matched by correlation in three consecutive keyframes. Fig. 6 shows the 187588 reconstructed points by incremental SfM for the complete sequence. The SfM drift is also visible thanks to an aerial photography (unused by our method). Fig. 7 shows the surface obtained at five different times t. The observer moves in the scene such that he/she is observing the most recent part of the surface at a (roughly) constant distance. This part is mainly within a ball whose center is c t-2 at t. At time t, the observer is located at c t-20 and is looking towards c t-2 . The observer and the surface end come forward simultaneously. These images are extracted from the joint video at http://maxime.lhuillier.free.fr.

Fig. 8 shows one global view and two local views of the last surface, which has 234354 triangles and 117145 vertices. 89.3% of free-space tetrahedra are outside tetrahedra. Remember that the surface is closed, so it also models the sky. In the figure, we remove the triangles on the sky to make viewing easier. Also we note that enforcing the man- ifold constraint is not an option. If we simply define the final surface S by the list of triangles between free-space and matter tetrahedra, we find that 25.5% of the S vertices are not regular. Fig. 1 shows that this degrades the quality. In the other cases, "Manifold" is between 50 and 600 ms. Thanks to Fig. 9 (on the right), we see that the computation times of "Manifold" and "Post-Processing"globally increase if td t increases and td t > 50. . Furthermore, td t < 280 in the whole sequence. Remember that d t is the smallest date of all outside tetrahedra destroyed at time t by "Delaunay". These results are consistent with those of a theoretical time complexity study: Delaunay and Carving are O(1), Manifold is O((td t ) log(td t )), Post-Processing is O(td t ). Actually, these tight bounds should be considered as conjectures since the proofs use strong assumptions and will be submitted in another paper.

We now explain the large values of "Manifold" if t ∈ [925, 1032]. In a complete trajectory loop, vertices added at the loop end (at time t) destroy outside tetrahedra created at the loop beginning (at time d t ) since these vertices and tetrahedra have similar 3d locations. The larger the loop, the larger td t , and the larger the "Manifold" (and "Post-Processing") computation time. Fig. 6 shows that the reconstructed trajectory has two incomplete (about 75%) loops: a large one on the top and a small one on the bottom. Here the loops are incomplete but the same principle applies for the small loop which is 75% closed if t ∈ [925, 1032]: there are times in [925, 1032] such that added vertices destroy outside tetrahedra created at the loop beginning. This does not apply in the large loop case since (1) the added vertices and outside tetrahedra are in a tubular neighborhood of the camera trajectory and (2) the neighborhood radius is less than the (divided by 2) distance between both ends of the loop. Fig. 8 shows the neighborhood and its size; the small loop is on the top and the large loop on the bottom.

Conclusion

To our knowledge, this paper presents the first system with four features: incremental reconstruction for triangulated manifold surface from sparse point cloud generated by SfM. In experiments, we use a synthetic image sequence to compare and discuss the performance of our incremental method and the related batch method. Although the batch method has slightly better results than our incremental method, the latter is interesting for online applications that the former can not solve. We also reconstruct parts of a city using a hand-held omnidirectional camera and provide Several steps of the method can be improved and are subjects for future work. Image edges could be reconstructed and integrated in the Delaunay to improve the surface. The region-growing step should be accelerated in the case of a closed loop in the camera trajectory. Currently, the surface is denoised assuming that the point cloud is dense enough to estimate a discrete Laplacian. This improves the surface but it would be better to design a dedicated denoising method for sparse SfM point clouds. Lastly, we plan to use our method for online applications, larger data sets, and to initialize surface reconstruction methods which are more time expensive and more accurate.

Appendix: Region Growing Acceleration

Here is a note to accelerate the "2-Manifold Extraction" step. In the algorithm of Section 3, we first insert tetrahedron ∆ in O and then check if δO is 2-manifold. However, we can do this faster (as in [START_REF] Boissonnat | Geometric structures for three-dimensional shape representation[END_REF]) if we first check a condition on the neighborhood of ∆ and then add ∆ to O (if the condition is meet). Let f be the number of ∆ facets which are in δO. If f = 1, ∆ is added to O if and only if the vertex of ∆, which is not in the δO facet, does not have adjacent tetrahedron in O. If f = 0, ∆ is added to O if and only if the four vertices meet this same condition. If f = 2, ∆ is added to O if and only if the edge of ∆, which is not an edge of the two δO facets, does not have adjacent tetrahedron in O. If f = 3 or f = 4, ∆ is added to O. In our implementation, we greatly accelerate these computations by precalculating for each vertex the list of its adjacent tetrahedra.

Figure 1

 1 Figure 1. (A) is a denoised non-manifold, (B) is a non-denoised 2-manifold, (C) is a denoised 2-manifold, (D) is the texture of (C). In (A,B,C), the triangle normals are encoded by colors.

Figure 2 .

 2 Figure 2. In the 2d case, surfaces, tetrahedra and triangles are replaced by curves, triangles and edges, respectively. (A) curves to be reconstructed and points which sample the curves. (B) the Delaunay triangulation of the points (the bold-faced Delaunay edges approximate the curves). (C) points (circle) and camera locations (squares) linked by rays (segments). (D) the Delaunay triangles intersected by rays are free-space and are denoted by "F". The other triangles are matter. The Delaunay edges separating free-space and matter are bold-faced.

  // **** initialization of priority queue (heap) Q **** Q = ∅; if (O==∅) { // used by batch algo. let ∆ ∈ F 0 be such that r(∆) is maximum; Q ← Q ∪ {∆}; } else // used by topology extension and incremental algo. for each tetrahedron ∆ in Q 0 ∩ F 0 if (∆ / ∈ O and one of its 4 neighbors is in O) Q ← Q ∪ {∆}; // **** region growing of O **** while (Q!=∅) { pick from Q the ∆ which has the largest r(∆); if (∆ ∈ O) continue; O ← O ∪ {∆}; if (all vertices of ∆ are regular) { // read the Appendix for each ∆ ′ in the list of the four ∆ neighbors

Figure 4 .

 4 Figure 4. Region growing for image t = 98. The number of tetrahedra layers contained in a pack is n = 20. Left: before region growing for image 98, the 2-manifolds of layers 20, 40, 60, 80 and 97 are already computed. Left-middle: point insertion destroys tetrahedra. The earliest creation date of tetrahedra which are destroyed due to point additions in T is dt = 51. Right-middle: the 2-manifolds of layers 60, 80, and 97 are invalid and destroyed. Right: region growing from layer ni0 to layer 98 by pack of n.

Figure 5 .

 5 Figure 5. Two omnidirectional images, the sparse point cloud (and camera trajectory) by SfM, the batch surface, the final incremental surface. We remove the triangles on the sky to make viewing in the figure easier. timated locations and the ground truth locations of the camera, respectively (c i is at the camera center and c g i is at the mirror apex). We found E(R)/600 = 5.1 cm and use R to map the estimated geometry (poses and point cloud) in the ground truth coordinate system. Both batch and incremental methods select points using ǫ = 10 degrees and add 2 extra points in the neighborhood of every c t (3d Delaunay Triangulation step). The incremental method has parameters k = 40 (Ray-Tracing step), n = 60 and b 0 = 10 (2-Manifold Extraction step), b 1 = 10 (Topology Extension step).

Figure 6 .

 6 Figure 6. From left to right: our hand-held camera, two images of the sequence, aerial view of the trajectory, the sparse point cloud reconstructed by incremental SfM.

Figure 7 .

 7 Figure 7. Images of the incremental surface reconstruction (also in the joint video). Top: gray levels encode the triangle normals. Bottom: one omnidirectional image is used for texture mapping. The black areas are due to triangles without texture in this image.

Figure 8 .

 8 Figure 8. Views of the incremental (final) surface.

Fig. 9 (

 9 on the left) shows the computation times of the different steps at each time t: "Delaunay" in yellow (3d Delaunay Triangulation+Dating), "Carving" in blue (Ray Tracing), "Manifold" in red (2-Manifold Ex-traction+Topology Extension), "Post-processing" in green (Surface Denoising) and "Total" in black. We use a Core 2 Duo E8500 at 3.16 GHz. About 121 points per time are added to the 3d Delaunay triangulation. "Delaunay" and "Post-processing" have almost negligible computation times in comparison to the other steps. "Carving" is less than 190 ms. If t ∈ [0, 925], "Manifold" is less than 200 ms.

Figure 9 .

 9 Figure 9. Calculation times of the incremental surface reconstruction as a function of t (left) or a function of tdt (right). a detailed explanation of the computation times.Several steps of the method can be improved and are subjects for future work. Image edges could be reconstructed and integrated in the Delaunay to improve the surface. The region-growing step should be accelerated in the case of a closed loop in the camera trajectory. Currently, the surface is denoised assuming that the point cloud is dense enough to estimate a discrete Laplacian. This improves the surface but it would be better to design a dedicated denoising method for sparse SfM point clouds. Lastly, we plan to use our method for online applications, larger data sets, and to initialize surface reconstruction methods which are more time expensive and more accurate.


Vadim Litvinov

Shuda Yu

Maxime Lhuillier

2-MANIFOLD RECONSTRUCTION FROM SPARSE VISUAL FEATURES

Keywords: 2-Manifold, Surface Reconstruction, 3D Delaunay Triangulation, Structure-from-Motion, Curves

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The automatic surface reconstruction of a scene from an image sequence is still an active topic. First, we summarize the advantages of the previous work [START_REF] Yu | Surface Reconstruction of Scenes using a Catadioptric Camera[END_REF] which estimates a 2-manifold surface directly from the sparse cloud of 3D points reconstructed by Structure-from-Motion (SfM). Then, we present our contributions.

The input point cloud is sparse (not dense) since only interest points, such as Harris or SIFT, are reconstructed. This greatly reduces the complexity of a surface reconstruction from points, and this is interesting for the compact 3D modeling of large scenes such as cities. Furthermore, a surface estimated from a sparse SfM cloud (and its visibility information) can be used as the initialization for a dense stereo step [START_REF] Hiep | Towards High-Resolution Large-Scale Multi-View Stereo[END_REF], which deforms the surface such that the global photo-consistency is optimized (this improves the surface quality). The dense stereo step is outside the paper scope.

The output surface is a 2-manifold; this is in contrast to other surface reconstruction methods [START_REF] Hilton | Scene Modeling from Sparse 3D Data[END_REF][START_REF] Labatut | Efficient Multi-View Reconstruction of Large-Scale Scenes using Interest Points, Delaunay Triangulation and Graph-Cut[END_REF][START_REF] Pan | ProFORMA: Probabilistic Feature-Based On-Line Rapid Model Acquisition[END_REF][START_REF] Lovi | Incremental Free-Space Carving for Real-Time 3D-Reconstruction[END_REF] with the same sparse input. For memory, a surface is a 2-manifold if and only if every point in the surface has a surface neighborhood which is topologically a disk. Since our surface is triangulated, every triangle is exactly connected by its 3 edges to 3 different triangles, the surface has neither holes nor selfintersections, and it cuts the 3D space into inside and outside regions. The manifold property is used by a lot of computer Thanks to CNRS, UBP, and Region of Auvergne for funding.

graphic methods [START_REF] Botsch | Polygon Mesh Processing[END_REF] and dense stereo methods which enforce smoothness or curvature constraints on the surface (e.g. [START_REF] Esteban | Silhouette and Stereo Fusion 3D Object Modeling[END_REF]).

The contributions of the paper are the following: we add reconstructed curves to the input point cloud in [START_REF] Yu | Surface Reconstruction of Scenes using a Catadioptric Camera[END_REF], and we experiment on both multiview image data with ground truth [START_REF] Seitz | A Comparison and Evaluation of Multi-View Stereo Reconstruction Algorithms[END_REF][START_REF] Strecha | On Benchmarking Camera Calibration and Multi-View Stereo for High Resolution Imagery[END_REF] and omnidirectional video image sequences (Ladybug). The work [START_REF] Faugeras | Representing Stereo Data with the Delaunay Triangulation[END_REF][START_REF] Morris | Image Consistent Surface Triangulation[END_REF][START_REF] Taylor | Surface Reconstruction from Feature-Based Stereo[END_REF][START_REF] Yu | Surface Reconstruction of Scenes using a Catadioptric Camera[END_REF] (including ours [START_REF] Yu | Surface Reconstruction of Scenes using a Catadioptric Camera[END_REF]) on 2-manifold reconstruction from sparse SfM data does not do that. Our curve reconstruction method involves matching based on correlation and assume that the epipolar constraint is known (this is different from curve reconstruction in [START_REF] Fabbri | 3D Curve Sketch: Flexible Curve-Based Stereo Reconstruction and Calibration[END_REF] where the correlation is not used). Section 2 summarizes the surface reconstruction method and describes the curve reconstruction. Sections 3 and 4 are the experiments and the conclusion.

AUTOMATIC SCENE MODELING

Here we briefly summarize two steps: SfM (Section 2.1) and surface reconstruction from points and visibilities (Section 2.2). Then Section 2.3 describes the curve reconstruction.

Sparse SfM

Our SfM is similar to [START_REF] Mouragnon | Generic and Real-Time Structure-From-Motion using Local Bundle Adjustment[END_REF] in the central generic case (the camera is calibrated) and alternates the steps below. First, a new key-frame is selected from the input video and interest points are matched with the previous key-frame using correlation. The new key-frame is such that the number of its matches with the two previous key-frames is larger than a threshold. Then, the new pose is robustly estimated (using Grunert's method [START_REF] Haralick | Review and analysis of solutions of the three point perspective pose estimation problem[END_REF] and RANSAC) and new 3D points are reconstructed from the new matches. Lastly, local bundle adjustment refines the geometry of the 3-most recent key-frames. The geometry includes both points and camera poses.

Surface Reconstruction from Points and Visibilities

This is the method in [START_REF] Yu | Surface Reconstruction of Scenes using a Catadioptric Camera[END_REF]. First, the 3D Delaunay triangulation T of the list P of reconstructed points is build: T partitions the convex hull of P , and the circumsphere (sphere circumscribing tetrahedron) of every tetrahedron in T does not contain any vertex of T within it. Second, every tetrahedron in T is labeled either free-space or matter with the help of the rays provided by SfM. A ray is a line segment connecting a 3D point to a viewpoint used to reconstruct it. A tetrahedron is free-space if it is intersected by a ray, otherwise it is matter. Third, the labeling is used to generate the surface. Several methods [START_REF] Pan | ProFORMA: Probabilistic Feature-Based On-Line Rapid Model Acquisition[END_REF][START_REF] Lovi | Incremental Free-Space Carving for Real-Time 3D-Reconstruction[END_REF][START_REF] Labatut | Efficient Multi-View Reconstruction of Large-Scale Scenes using Interest Points, Delaunay Triangulation and Graph-Cut[END_REF] directly consider the target surface as the list of triangles separating the free-space and matter. Unfortunately, the resulting surface can be non-manifold.

In order to generate a 2-manifold, a region-growing approach is applied. The idea is to grow iteratively an outside region O in the free-space region by maintaining the manifold property on the border of O during the growing [START_REF] Yu | Surface Reconstruction of Scenes using a Catadioptric Camera[END_REF]. Before the growing, we assign inside to all tetrahedrons. The final manifold is the list of triangles separating the outside tetrahedrons and inside tetrahedrons.

Lastly, post-processing is used to improve the quality of the surface. It includes peak removal, surface denoising using discrete Laplacian, triangle removal in the sky (optional), and surface texturing.

Reconstruct and Integrate Curves

Here we reconstruct curves in 3D from image edges. The result is a list of sampled 3D points along the curves, which are added to the SfM points as inputs to the surface estimation. There are several steps: detection and matching in Section 2.3.1, reconstruction and sampling in Section 2.3.2.

Detection and Matching

First of all, the edge curves are detected in the incoming keyframe using the well known Canny operator. Then, they are parametrized in the horizontal or vertical direction (depending on the main curve orientation), refined using the gradient orientation information provided by Canny, and filtered to remove the curves that are too short for being useful.

Then the curves are matched with those in the previous keyframe. Unfortunately, curve geometry based matching method as described in [START_REF] Fabbri | 3D Curve Sketch: Flexible Curve-Based Stereo Reconstruction and Calibration[END_REF] exhibited bad performance when the images contained a lot of condensed and similar curves (for example, as in the Ladybug datasets, see Section 3), it just can't make a matching decision based on curve geometry alone. So, we decided to use a correlation based method heavily inspired by a classical points matching.

To perform a curve matching between images I 1 and I 2 we proceed as follow: for each curve C 1 in I 1 we construct an epipolar lines strip in the image I 2 , then we reject all the curves not at least partially included into it. The remaining curves of I 2 form a list of candidates L c . For each curve C 2 in L c we compute a zero-normalized cross-correlations (ZNCC) between intersections of C 2 with the epipolar lines and the points of C 1 that produced them. Finally, the mean value of all these correlations define the score of the curve C 2 . C 1 is matched with the curve of L c having the best score.

Reconstruction and Sampling

First, the lists of curve tracks generated in Section 2.3.1 are converted to lists of image point tracks thanks to the epipolar constraints. A curve track is a list of matched curves c i0 , c i0+1 , ...c i1 such that c i is in the i-th key-frame. A point track is a list of matched points p i0 , p i0+1 , ...p i1 where p i is in the i-th key-frame: p i+1 is obtained from p i by the intersection of c i+1 and the epipolar line of p i in the i + 1-image. If the angle between the epipolar line and the c i+1 tangent at p i+1 is small, the intersection is inaccurate and the point track is truncated (a point track can be shorter than the curve track).

Then, every point track is reconstructed by ray intersection: RANSAC is used for robustness and initializing a Levenberg-Marquardt optimization [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. Unfortunately, using this direct approach leads to some strange artifacts (left of Fig. 1 for "Temple" [START_REF] Seitz | A Comparison and Evaluation of Multi-View Stereo Reconstruction Algorithms[END_REF] dataset). In fact, all the points of the same curve do not necessarily use the same views to initialize themselves by RANSAC. This difference in initial position leads to zigzag like 3D curves. To avoid this problem, we perform a two-step initialization: first, all the points are initialized independently; then, at the second stage, we select the two views that were used the most often and all the points are initialized using these two views. As can be seen on the right image of Fig. 1, this approach greatly improves the reconstructed curves.

Finally, we only retain for surface estimation one reconstructed point over s points for every curve track. The step size s is large enough to preserve speed and sparsity of the global algorithm and at the same time small enough to ensure that curves are correctly integrated into the Delaunay triangulation. (c) Ground truth.

Fig. 2. "Temple" dataset reconstruction results.

EXPERIMENTS

Our method is evaluated on three different datasets: "Temple" dataset [START_REF] Seitz | A Comparison and Evaluation of Multi-View Stereo Reconstruction Algorithms[END_REF], "Fountain-P11" and "Herzjesu-P8" dataset [START_REF] Strecha | On Benchmarking Camera Calibration and Multi-View Stereo for High Resolution Imagery[END_REF], and finally our own "Hall" and "Aubière" dataset. The first two are provided with the ground truth so this allows us to quantitatively evaluate the performance of our algorithm and to compare it with the others. The third allows us to qualitatively evaluate our method on a cluttered interior and a residential scene. All the processing times are evaluated on a 4xIntel Xeon W3530 at 2.8 GHz. A summary of important values about every dataset presented in this section can be found in table 1. Note that our method was initially designed for the reconstruction of complete environments using omnidirectional cameras as in our own (third) dataset. In this case, the camera positions are assumed to be inside the convex hull of the reconstructed scene points. In the other cases (first and second datasets) a workaround is needed: we add Steiner vertices (extra points) at the border of a large bounding box of the object ("Temple") or behind the camera ("Fountain-Herzjesu") to the Delaunay triangulation, then we apply our method and remove the triangles which are incident to these Steiner vertices.

"Temple" dataset

The "Temple" provided by [START_REF] Seitz | A Comparison and Evaluation of Multi-View Stereo Reconstruction Algorithms[END_REF] is a standard multiview reconstruction dataset (312 separate views). In this evaluation, we use the camera positions provided with the dataset to reconstruct the 3D points. Two models were sent for evaluation to the dataset creators. The first (see Fig. 2.a) is reconstructed using the points only and the second one (see Fig. This lead us to two conclusions. Firstly, the use of the curves provides indeed an enhancement to the resulting model. Secondly, compared to other methods evaluated using the same datasets, our algorithm is, from the best of our knowledge, the fastest CPU based reconstruction method, and the resulting precision is not so bad compared to majority of the dense stereo algorithms. This makes our method a good candidate for initialization of a dense stereo method.

"Fountain-P11" and "Herzjesu-P8" datasets

"Fountain" and "Herzjesu" are another standard set of multiview stereo datasets provided with their respective ground truths (courtesy of [START_REF] Strecha | On Benchmarking Camera Calibration and Multi-View Stereo for High Resolution Imagery[END_REF]). They contain, respectively, 11 and 8 high-resolution images. In this evaluation, we use the camera positions provided with the dataset to reconstruct the 3D points.

To evaluate the error distribution against the ground truth we subdivide our output mesh with a very thin step and, for each points of the output mesh, we calculate the distance to the nearest point of the ground truth mesh.

The results can be seen on Figs. 3 and4 respectively. Part (a) provides an error distribution cumulative histogram. Part (b) shows the output mesh and part (c) displays the estimated model colored in function of the error value (green is better, scaled between values at 20% and 80% for better visibility). The computations of "Fountain" and "Herzjesu" take 123 s. and 74 s., respectively.

As can be seen, our method provides honorable results, so it can be a good choice when the speed of computation is at least as important as precision. We do not show the results obtained by mixing points with curves because adding curves provides a very minor improvement. "Hall" and "Aubière" datasets are two large datasets taken in and in the proximity of Institut Pascal with a PointGrey Ladybug omnidirectional camera. Ladybug is a rigid multicamera system consisting of six synchronized cameras each of which takes 1024x768 images at 15 frames/second. Unfortunately, ground truth is not available for these acquisitions, so we provide only qualitative results. This allows us to evaluate the general performance of the complete algorithm as described in the Section 2. Fig. 5 shows the complete Ladybug frames for both sequences.

"Hall" sequence is roughly 20 m. long and contains 1211 frames (6 images each). There are 117 selected key-frames; the total processing time is about 19 min. including curves (the main computation part is due to SfM). Reconstruction of the representative scene can be seen on Fig. 6: part (b) contains a real image of the scene for reference. We see that mixing points and curves provides a better result than the usage of the points alone.

"Aubière" sequence is around 700 m. long and contains 3140 frames. There are 655 selected key-frames and the total processing time is about 56 min. (the main computation part is due to SfM). The results are on Fig. 7. Contrary to the previous sequence, the curves doesn't bring any significant improvement over the points alone. Global views without curves are also shown in Fig. 8;

CONCLUSION

The contributions of this paper are (1) the evaluation of the 2-manifold surface sparse reconstruction method described in [START_REF] Yu | Surface Reconstruction of Scenes using a Catadioptric Camera[END_REF] against ground truth in some common multiview stereo datasets and (2) the evaluation of the enhancement provided by mixing the input points cloud of this algorithm with reconstructed curves.

According to the comparison of the output surface with the ground truth, the precision achieved by this method cannot yet compete with the dense stereo methods, but we think that it is sufficient to be used as initialization of dense stereo. Furthermore, given the achieved speed, the algorithm is indeed a good choice not only for this purpose but also for the applications where the speed is at least as important as the precision.

On the other hand, mixing the curves with the points of interest provides some significant enhancement if input images are only slightly textured as in "Temple" and "Hall". But it is of little interest in the textured case, as for "Fountain" and "Herzjesu". So there are no general conclusion and the real benefits must be evaluated for each sequence individually. Nevertheless, it will be interesting to perform more investigations in this matter as better results can certainly be achieved.

  (a) Classical RANSAC. (b) Unified RANSAC.

Fig. 1 .

 1 Fig. 1. Different 3D curve initialization methods.

  (a) Points only. (b) Points and curves.

  2.b) is reconstructed using points and curves. The ground truth image is provided for reference on Fig. 2.c. The point-only model is reconstructed in 31 s. and provides an accuracy of 0.66 mm. (for 90% of reconstructed points) and completeness of 93% (for 1.25 mm. error). The points and curves reconstruction takes 59 s. and provides an accuracy of 0.59 mm. and completeness of 95.4%. (These values can be compared to those of other algorithms at the following address:http://vision.middlebury.edu/ mview/eval/).

  of total) (a) Error distribution. (b) Reconstructed surface. (c) Ground truth surface with error shown in pseudocolor.

Fig. 3 .

 3 Fig. 3. "Fountain" dataset reconstruction results.

  of total) (a) Error distribution. (b) Reconstructed surface. (c) Ground truth surface with error shown in pseudocolor.

Fig. 4 .

 4 Fig. 4. "Herzjesu" dataset reconstruction results.

3. 3 .

 3 "Hall" and "Aubière" datasets (a) "Hall". (b) "Aubière".

Fig. 5 .

 5 Fig. 5. Ladybug images for "Hall" and "Aubière" video sequences.

  (a) SfM point cloud and Ladybug trajectory. (b) Viewed part of the "Hall". (c) Scene reconstructed using points only. (d) Scene reconstructed using points and curves.

Fig. 6 .

 6 Fig. 6. "Hall" dataset reconstruction results.

  (a) SfM point cloud. (b) Textured version of the resulting mesh using points only. (c) Resulting mesh using points only. (d) Resulting mesh using points and curves.

Fig. 7 .

 7 Fig. 7. "Aubière" dataset reconstruction results.

Fig. 8 .

 8 Fig. 8. Global views of the "Aubière" surface using points only (best viewed in colors).

Table 1 .

 1 Descriptions of datasets used in our experimental tests (1k=1000).

		Temple	Fountain-P11 Herzjesu-P8	Hall	Aubière
	Image resolution	1x640x480	1x3072x2048 1x3072x2048	6x1024x768	6x1024x768
	Key-frames (total frames)	312 (312)	11 (11)	8 (8)	117 (1211)	493 (2443)
	Harris points per frame	1.5k	60k	50k	24k	24k
	curve sampling size s	4	8	8	8	8
	SfM 3D Points (curves)	32k (56k)	67k	29k	86k (145k)	295k (374k)
	Final triangles (curves)	49k (82k)	121k	52k	83k (305k)	399k (520k)
	Time: 3D points (+curves) 13 s. (+28 s.)	100 s.	64 s.	17 min. (+102 s.) 55 min. (+8 min.)
	Time: Surface	18 s.	23 s.	10 s.	47 s.	81 s.