Stéphane Zampelli

Yves Deville

Christine Solnon
email: christine.solnon@liris.cnrs.fr

Sébastien Sorlin
email: sebastien.sorlin@liris.cnrs.fr

Pierre Dupont

Filtering for Subgraph Isomorphism

A subgraph isomorphism problem consists in deciding if there exists a copy of a pattern graph in a target graph. We introduce in this paper a filtering algorithm dedicated to this problem. The main idea is to label every node with respect to its relationships with other nodes of the graph, and to define a partial order on these labels in order to express compatibility of labels for subgraph isomorphism. This partial order over labels is used to filter domains. Labelings can also be strengthened by adding information from the labels of the neighbors. Such a strengthening can be applied iteratively until a fixpoint is reached. Practical experiments illustrate that our new filtering approach is more effective on difficult instances of scale free graphs than state-of-the-art algorithms and other CP approaches.

Introduction

Graphs are widely used in real-life applications to represent structured objects, e.g., molecules, images, or biological networks. In many of these applications, one looks for a copy of a pattern graph into a target graph [START_REF] Conte | Thirty years of graph matching in pattern recognition[END_REF]. This problem, known as subgraph isomorphism, is NP-complete [START_REF] Garey | Computers and Intractability[END_REF] in the general case.

There exists dedicated algorithms for solving subgraph isomorphism problems, such as [START_REF] Ullmann | An algorithm for subgraph isomorphism[END_REF][START_REF] Cordella | An improved algorithm for matching large graphs[END_REF]. However, such dedicated algorithms can hardly be used to solve more general problems, with additional constraints, or approximate subgraph isomorphism problems, such as the one introduced in [START_REF] Zampelli | Approximate constrained subgraph matching[END_REF].

An attractive alternative to these dedicated algorithms is Constraint Programming (CP), which provides a generic framework for solving constraint satisfaction problems (CSP). Indeed, subgraph isomorphism problems may be formulated as CSP in a straightforward way [START_REF] Régin | Développement d'Outils Algorithmiques pour l'Intelligence Artificielle[END_REF][START_REF] Rudolf | Utilizing constraint satisfaction techniques for efficient graph pattern matching[END_REF]. To make CP competitive with dedicated approaches for these problems, [START_REF] Larrosa | Constraint satisfaction algorithms for graph pattern matching[END_REF] has introduced a global monomorphism constraint, and an associated filtering algorithm, together with redundant Alldiff constraints. [START_REF] Zampelli | Approximate constrained subgraph matching[END_REF] has extended this work to approximate subgraph isomorphism, and has shown that CP is competitive with dedicated approaches.

Contribution. In this paper, we introduce a new filtering algorithm for the subgraph isomorphism problem that exploits the global structure of the graph to achieve a stronger partial consistency. This work takes inspiration from the partition refinement procedure used in Nauty [START_REF] Mckay | Practical graph isomorphism[END_REF] and Saucy [START_REF] Darga | Exploiting structure in symmetry detection for cnf[END_REF] for finding graph automorphisms: the idea is to label every node by some invariant property, such as node degrees, and to iteratively extend labels by considering labels of adjacent nodes. Similar labelings are used in [START_REF] Sorlin | A global constraint for graph isomorphism problems[END_REF][START_REF] Sorlin | A new filtering algorithm for the graph isomorphism problem[END_REF] to define filtering algorithms for the graph isomorphism problem: the idea is to remove from the domain of a variable associated to a node v every node the label of which is different from the label of v. The extension of such a label-based filtering to subgraph isomorphism problems mainly requires to define a partial order on labels in order to express compatibility of labels for subgraph isomorphism: this partial order is used to remove from the domain of a variable associated to a node v every node the label of which is not compatible with the label of v. We show that this extension is more effective on difficult instances of scale free graphs than state-of-the-art subgraph isomorphism algorithms and other CP approaches.

Outline of the paper. Section 2 presents the subgraph isomorphism problem and related CP models. Section 3 describes the theoretical framework of our filtering: it first introduces the concept of labeling, and shows how labelings can be used for filtering; it then shows that labelings can be iteratively strengthened by adding information from labels of neighbors. Section 4 introduces the practical framework and describes how to compute a label strengthening. An exact algorithm as well as an approximate version are provided. Experimental results are described in Section 5.

Subgraph Isomorphism

Definitions

An (undirected) graph G = (N, E) consists of a node set N and an edge set E ⊆ N × N , where an edge (u, v) is a couple of nodes.

A subgraph isomorphism problem between a pattern graph G p = (N p , E p ) and a target graph G t = (N t , E t ) consists in deciding whether G p is isomorphic to some subgraph of G t . More precisely, one should find an injective function

f : N p → N t such that ∀(u, v) ∈ N p × N p , (u, v) ∈ E p ⇒ (f (u), f (v)) ∈ E t .
The problem is also called subgraph monomorphism problem or subgraph matching in the literature. The function f is called a subgraph matching function.

In the following, we assume G p = (N p , E p ) and G t = (N t , E t ) to be the underlying instance of subgraph isomorphism problem. We also define N ode = N p ∪ N t , Edge = E p ∪ E t , n p = #N p , n t = #N t , n = #N ode, d p and d t the maximal degree of the graphs G p and G t , and d = max(d p , d t ).

CP Models for Subgraph Isomorphism

A subgraph isomorphism problem can be formulated as a CSP in a straightforward way [START_REF] Régin | Développement d'Outils Algorithmiques pour l'Intelligence Artificielle[END_REF][START_REF] Rudolf | Utilizing constraint satisfaction techniques for efficient graph pattern matching[END_REF][START_REF] Larrosa | Constraint satisfaction algorithms for graph pattern matching[END_REF]. A variable x u is associated with every node u of the pattern graph and its domain is the set of target nodes. A global Alldiff constraint [START_REF] Regin | A filtering algorithm for constraints of difference in CSPs[END_REF] ensures that the matching function is injective. Edge matching is ensured by a set of binary constraints:

∀(u, v) ∈ N p × N p , c2(x u , x v ) ≡ ((u, v) ∈ E p ⇒ (x u , x v ) ∈ E t ) .

Theoretical Framework

This section introduces a new filtering algorithm for subgraph isomorphism. We will show in the next section how filtering can be achieved in practice from this theoretical work.

Subgraph Isomorphism Consistent Labelings

Definition 1. A labeling l is defined by a triple (L, , α) such that -L is a set of labels that may be associated to nodes; -⊆ L × L is a partial order on L; α : Node → L is a total function assigning a label α(v) to every node v.

A labeling induces a compatibility relation between nodes of the pattern graph and the target graph. Definition 2. The set of compatible couples of nodes induced by a labeling l = (L, , α) is defined by

CC l = {(u, v) ∈ N p × N t | α(u) α(v)}
This compatibility relation can be used to filter the domain of a variable x u associated with a node u of the pattern graph by removing from it every node v of the target graph such that (u, v) ∈ CC l .

The goal of this work is to find a labeling that filters domains as strongly as possible without removing solutions to the subgraph isomorphism problem, i.e., if a node v of the pattern graph may be matched to a node u of the target graph by a subgraph matching function, then the label of v must be compatible with the label of u. This property is called subgraph isomorphism consistency.

Definition 3. A labeling l is subgraph isomorphism consistent (SIC) iff for any subgraph matching function f , we have ∀v ∈ N p , (v, f (v)) ∈ CC l .
In the context of graph isomorphism, such as in [START_REF] Mckay | Practical graph isomorphism[END_REF], as opposed to subgraph isomorphism studied here, an SIC labeling is often called an invariant. In this case, the partial ordering is replaced by an equality condition: two nodes are compatible if they have the same label.

Many graph properties, that are "invariant" to subgraph isomorphism, may be used to define SIC labelings such as, e.g., the three following SIC labelings:

l deg = (N, ≤, deg) where deg is the function that returns node degree; l distance k = (N, ≤, distance k ) where distance k is the function that returns the number of nodes that are reachable by a path of length smaller than k; l clique k = (N, ≤, clique k ) where clique k is the function that returns the number of cliques of size k that contains the node. Example. Let us consider for example the subgraph isomorphism problem displayed in Fig. 1. Note that this instance has no solution as G p cannot be mapped into a subgraph of G t . The labeling l deg = (N, ≤, deg) assigns the following labels to nodes.

deg(A) = deg(B) = deg(D) = deg(2) = deg(4) = 4 deg(C) = deg(E) = deg(F ) = deg(G) = deg(1) = deg(3) = 3 deg(5) = deg(6) = 2
Hence, the set of compatible couples induced by this labeling is

CC l deg = {(u, v) | u ∈ {2, 4}, v ∈ {A, B, D}} ∪ {(u, v) | u ∈ {1, 3, 5, 6}, v ∈ {A, B, C, D, E, F, G}}
This set of compatible couples allows one to remove values C, E, F and G from the domains of the variables associated with nodes 2 and 4.

Strengthening a Labeling

We propose to start from an elementary SIC labeling that is easy to compute such as the l deg labeling defined above, and to iteratively strengthen this labeling. The strength of labelings is defined with respect to the induced compatible couples as follows.

Definition 4. Let l and l be two labelings. l is strictly stronger than l iff CC l ⊂ CC l , and l is equivalent to l iff CC l = CC l .

A stronger labeling yields a better filtering, as it contains less compatible couples.

To strengthen a labeling, the idea is to extend the label of a node by adding information from the labels of its neighbors. This information is a multiset (as several neighbors may have the same label). We shall use the following notations for multisets. Definition 5. Given an underlying set A, a multiset is a function m : A → N, such that m(a) is the multiplicity (i.e., the number of occurrences) of a in m. The multiset m can also be represented by the bag {a 0 , . . . , a 0 , a 1 , . . .} where elements are repeated according to their multiplicity.

For example, the multiset m that contains 2 occurrences of a, 3 occurrences of b, and 1 occurrence of c is defined by m(a)=2, m(b)=3, m(c)=1, and ∀x ∈ {a, b, c}, m(x)=0. This multiset may also be represented by {a, a, b, b, b, c}.

Given a partial order on a set A, we extend the partial order to multisets over A as follows. Definition 6. Given two multisets m and m over a set A, and a partial order ⊆ A×A, we define m m iff there exists a total injective mapping t : m → m such that ∀a i ∈ m, a i t(a i ).

In other words, m m iff for every element of m there exists a different element of m which is greater or equal. For example, if we consider the classical ordering on N, we have {3, 3, 4} {2, 3, 5, 5}, but {3, 3, 4} is not comparable with {2, 5, 6}. Note that comparing two multisets is not trivial in the general case, especially if the order relation on the underlying set A is not total. This point will be handled in the next section.

We now define the labeling extension procedure.

Definition 7. Given a labeling l = (L, , α), the neighborhood extension of l is the labeling l = (L , , α ) such that:

every label of L is composed of a label of L and a multiset of labels of L, i.e., L = L • (L → N); the labeling function α extends every label α(v) by the multiset of the labels of the neighbors of v, i.e., α

(v) = α(v) • m where ∀l i ∈ L, m(l i ) = #{u | (u, v) ∈ Edge ∧ α(u) = l i }; -the partial order on the extended labels of L is defined by l 1 • m 1 l 2 • m 2 iff l 1 l 2 and m 1 m 2 .
The next theorem states that the neighborhood extension of a SIC labeling is a stronger (or equal) SIC labeling.

Theorem 1. Let l = (L, , α) be a labeling, and l = (L , , α ) be its neighborhood extension. If l is an SIC labeling, then (i) l is also SIC, and (ii) l is stronger than or equal to l.

Proof. (i): Let f be a subgraph matching function and v ∈ N p . We show that α

(v) α (f (v)), that is α(v) α(f (v)) and m m , with m (resp. m ) the multiset of the labels of the neighbors of v in G p (resp. of f (v) in G t ): -α(v) α(f (v)) because l is SIC; -m m because m contains, for each neighbor u of v, the label α(f (u)) of
the node matched to u by the subgraph matching function f ; as l is SIC, α(u) α(f (u)), and thus m m .

(ii) : This is a direct consequence of the partial order on the extended labels in

L (Definition 7) : α(u) α(v) is one of the conditions to have α (u) α (v).
Example. Let us consider again the subgraph isomorphism problem displayed in Fig. 1, and the labeling l deg = (N, ≤, deg) defined in 3.1. The neighborhood extension of l deg is the labeling l = (L , , α ) displayed below. Note that we only display compatibility relationships l i l j such that l i is the label of a node of the pattern graph and l j is the label of a node of the target graph as other relations are useless for filtering purposes.

α (A) = 4 • {3, 3, 4, 4} α (B) = α (D) = 4 • {3, 3, 3, 4} α (2) = α (4) = 4 • {2, 2, 3, 3} 4 • {3, 3, 4, 4} and 4 • {3, 3, 3, 4} α (C) = 3 • {4, 4, 4} α (E) = α (F ) = 3 • {3, 4, 4} 4 • {3, 3, 4, 4}, 3 • {4, 4, 4} and 3 • {3, 4, 4} α (1) = α (3) = 3 • {3, 4, 4} 4 • {3, 3, 4, 4}, 3 • {4, 4, 4} and 3 • {3, 4, 4} α (G) = 3 • {3, 3, 4} α (5) = α (6) = 2 • {4, 4} 4 • {3, 3, 4, 4}, 3 • {4, 4, 4} and 3 • {3, 4, 4}
Hence, the set of compatible couples induced by this extended labeling is

CC l = {(u, v) | u ∈ {2, 4}, v ∈ {A, B, D}} ∪ {(u, v) | u ∈ {1, 3, 5, 6}, v ∈ {A, C, E, F }}
As compared to the initial labeling l deg , this set of compatible couples allows one to further remove values B, D and G from the domains of the variables associated with nodes 1, 3, 5 and 6.

Iterative Labeling Strengthening

The strengthening of a labeling described in the previous section can be repeated by relabeling nodes iteratively, starting from a given SIC labeling l.

Definition 8. Let l = (L, , α) be an initial SIC labeling. We define the sequence of SIC labelings l i = (L i , i , α i ) such that l 0 = l and l i+1 = neighborhood extension of l i (i ≥ 0).

A theoretical filter can be built on this sequence. Starting from an initial SIC labeling function l = l 0 , we iteratively compute l i+1 from l i and filter domains with respect to the set of compatible couples induced by l i until either a domain becomes empty (thus indicating that the problem has no solution) or reaching some termination condition.

A termination condition is to stop iterating when the sequence reaches a fixpoint, i.e., a step where any further relabeling cannot change the strength of the labeling. We show in [START_REF] Zampelli | Filtering for subgraph matching[END_REF] that such a fixpoint is reached in O(n p • n t ) steps, when both the set of compatible couples and the number of different labels are not changed between two consecutive steps.

Example. Let us consider again the subgraph isomorphism problem displayed in Fig. 1, and let us suppose that the sequence of SIC labelings is started from l 0 = l deg = (N, ≤, deg) as defined in 3.1. After the first iteration, the neighborhood extension l 1 of l 0 is the labeling displayed in the example of section 3.2. To improve reading, we rename these labels as follows:

α 1 (A) = 4 • {3, 3, 4, 4} renamed m 1 α 1 (B) = α 1 (D) = 4 • {3, 3, 3, 4} renamed m 2 α 1 (2) = α 1 (4) = 4 • {2, 2, 3, 3} renamed m 3 1 {m 1 , m 2 } α 1 (C) = 3 • {4, 4, 4} renamed m 4 α 1 (E) = α 1 (F ) = α 1 (1) = α 1 (3) = 3 • {3, 4, 4} renamed m 5 1 {m 1 , m 4 } α 1 (G) = 3 • {3, 3, 4} renamed m 6 α 1 (5) = α 1 (6) = 2 • {4, 4} renamed m 7 1 {m 1 , m 4 , m 5 }
From labeling l 1 , we compute the following extended labels and partial order:

α 2 (A) = m 1 • {m 2 , m 2 , m 4 , m 5 } renamed n 1 α 2 (B) = m 2 • {m 1 , m 4 , m 5 , m 6 } renamed n 2 α 2 (C) = m 4 • {m 1 , m 2 , m 2 } renamed n 3 α 2 (D) = m 2 • {m 1 , m 4 , m 5 , m 5 } renamed n 4 α 2 (E) = m 5 • {m 1 , m 2 , m 6 } renamed n 5 α 2 (F ) = m 5 • {m 2 , m 2 , m 6 } renamed n 6 α 2 (G) = m 6 • {m 2 , m 5 , m 5 } renamed n 7 α 2 (1) = α 2 (3) = m 5 • {m 3 , m 3 , m 5 } renamed n 8 2 {n 1 , n 3 } α 2 (2) = α 2 (4) = m 3 • {m 5 , m 5 , m 7 , m 7 } renamed n 9 2 {n 4 } α 2 (5) = α 2 (6) = m 7 • {m 3 , m 3 } renamed n 10 2 {n 1 , n 3 , n 5 , n 6 }
The set of compatible couples induced by this labeling is

CC l 2 = {(2, D), (4, D), (1, A), (1, C), (3, A), (3, C)} ∪ {(u, v) | u ∈ {5, 6}, v ∈ {A, C, E, F }}
This set of compatible couples allows one to remove values A and B from the domains of the variables associated with nodes 2 and 4. Hence, the domains of these two variables only contain one value (D), and thanks to the alldiff constraint on the variables, an inconsistency is detected.

Algorithm 1: Filtering procedure Input: two graphs Gp = (Np, Ep) and Gt = (Nt, Et) such that Np ∩ Nt = ∅, an initial SIC labeling l 0 = (L 0 , α 0 , 0 ), initial domains D : Np → P(Nt), a limit k on the number of iterations Output: filtered domains for every node u ∈ Np do:

D(u) ← D(u) ∩ {v ∈ Nt|α 0 (u) 0 α 0 (v)} 1 i ← 1 2 while ∀u ∈ Np, D(u) = ∅ and i ≤ k and fixpoint not reached do 3 for every node u ∈ Np ∪ Nt do 4 m i u ← multiset containing an occurrence of α i-1 (v), ∀(u, v) ∈ Ep ∪ Et 5 α i (u) ← α i-1 (u) • m i u 6 L i ← {α i (u)|u ∈ Np ∪ Nt}; rename labels in L i and α i 7 i ← {(α i (u), α i (v)) | u ∈ Np ∧ v ∈ D(u) ∧ test(m i u , m i v , i-1 )} 8 for every node u ∈ Np do: D(u) ← D(u) ∩ {v ∈ Nt|α i (u) i α i (v)} 9 i ← i + 1 10 return D 11 4 Practical Framework
Algorithm 1 describes the overall filtering procedure. Starting from an initial SIC labeling, that may be, e.g., l deg , this procedure first filters domains with respect to this initial labeling (line 1) and then iteratively extends this labeling (lines 4-8) and filters domains with respect to the new extended labeling (line 9) until some domain becomes empty, or a maximum number of iterations have been performed, or a fixpoint is reached.

Labeling extension (lines 4-8) is decomposed into three steps:

lines 4-6: α i is computed from α i-1 ; this step is done in O(#Edge); line 7: labels of L i are renamed; this step is done in O(d • #N ode); line 8: the partial order i is computed, i.e., for every couple of nodes (u, v) such that u is a node of the pattern graph and v is a node of the target graph which was compatible with u at step i -1, we test for the compatibility of the multisets m i u and m i v to determine if the labels of u and v are still compatible at step i. Testing the compatibility of two mulisets is not trivial. We show in 4.1 how to do this exactly in O(d 5/2 ), so that line 8 has a time complexity of O(n p • n t • d 5/2 ). We then show in 4.2 how to compute an order inducing a weaker filtering in O(n t • d • (n p + d t • log n t )). These two variants are experimentally compared in Section 5.

The filtering step (line 9) is done in O(n p • n t ).

Exact computation of the partial order

Given two multisets m u and m v , and a partial order , the function test(m u , m v ,

) determines if m u m v , i.e., if there exists for each label occurrence in m u a distinct label occurrence in m v which is greater or equal according to .

Property 1. Let G = (N = (N u , N v ), E) be the bipartite graph such that N u (resp. N v ) associates a different node with every label occurrence in the multiset m u (resp. m v ), and E contains the set of edges (i, j) such that i j. We have m u m v iff there exists a matching that covers N u in G.

Hopcroft [START_REF] Hopcroft | An n 5/2 algorithm for maximum matchings in bipartite graphs[END_REF] proposes an algorithm for solving this problem in

O(|N u | • |N v | • |N |).
As the sizes of m u and m v are bounded by the maximal degree d, the test function can be done in O(d 5/2 ).

Computation of an approximated order

If is a total order, the function test(m u , m v , ) can be implemented more efficiently, by sorting each multiset and matching every label of m u with the smallest compatible label of m v . In this case, the complexity of test is O(d•log d).

When is not a total order, one may extend it into a total order ≤. This total order can then be used in the test function to determine if m u ≤ m v . However, the total order introduces new label compatibilities so that test(m u , m v , ≤) may return true while test(m u , m v , ) returns false. As a consequence, using this approximated order may induce a weaker filtering.

In this section, we first introduce the theoretical framework that defines a new neighborhood labeling extension based on a total order and proves its validity; then we show how it can be achieved in practice.

Neighborhood labeling extension based on a total order. The next definition gives a simple condition on the total order to ensure its consistency with respect to the partial order, i.e., to ensure that test(m u , m v , ) = True ⇒ test(m u , m v , ≤) = True. Definition 9. Let l = (L, , α) be a labeling. A consistent total order for l is a total order ≤ on L such that ∀u ∈ n p , ∀v ∈ n

t , α(u) α(v) ⇒ α(u) ≤ α(v)
We extend the order ≤ on multisets like for partial orders in Definition 6, i.e., m ≤ m iff there exists an injective function t : m → m such that ∀a i ∈ m, a i ≤ t(a i ). Hence, m m ⇒ m ≤ m . Let us note however that this extension of ≤ to multisets only induces a partial order on multisets as some multisets may not be comparable.

We can then define a new neighborhood extension procedure, based on a consistent total order. Definition 10. Let l = (L, , α) be a labeling, and ≤ be a consistent total order for l. The neighborhood extension of l based on ≤ is the labeling l ≤ = (L , ≤ , α ) where L and α are defined like in Definition 7, and the order relation

≤ ⊆ L × L is defined by l 1 • m 1 ≤ l 2 • m 2 iff l 1 l 2 ∧ m 1 ≤ m 2
The next theorem shows that the neighborhood extension l ≤ based on ≤ may be used in our iterative labeling process, and that it is stronger or equal to l. However, it may be weaker than the neighborhood extension based on the partial order . Indeed, the total order induces more compatible couples of labels than the partial order.

Theorem 2. Let l = (L, , α), l = (L , , α ), and l ≤ = (L , ≤ , α ), be three labelings such that l is the neighborhood extension of l and l ≤ is the neighborhood extension of l based on a consistent total order ≤.

If l is an SIC labeling, then (i) l ≤ is SIC, (ii) l ≤ is stronger than (or equal to) l, and (iii) l is stronger than (or equal to) l ≤ .

Proof. (ii) and (iii): For labeling l , we have

l 1 •m 1 l 2 •m 2 iff l 1 l 2 ∧m 1 m 2 . As ≤ is consistent w.r.t. , we have m m ⇒ m ≤ m . Hence, CC l ⊆ CC l ≤ ⊆ CC l . (i) is a direct consequence of (iii), as l is SIC (Theorem 1).
Different consistent total orders may be derived from a given partial order, leading to prunings of different strength: the less new couples of compatible nodes are introduced by the total order, the better the filtering. However, we conjecture in [START_REF] Zampelli | Filtering for subgraph matching[END_REF] that finding the best consistent total order is NP-hard. Hence, we propose a heuristic algorithm that aims at computing a total order that introduces few new compatible couples without guarantee of optimality. Let us note L p (resp. L t ) the set of labels associated with nodes of the pattern graph G p (resp. target graph G t ). We shall suppose without loss of generality1 that L p ∩ L t = ∅. The idea is to sequence the labels of L p ∪ L t , thus defining a total order on these labels, according to the following greedy principle: starting from an empty sequence, one iteratively adds some labels of L p ∪ L t at the end of the sequence, and removes these labels from L p and L t , until L p ∪ L t = ∅.

To choose the labels added in the sequence at each iteration, our heuristic is based on the fact that the new couples of compatible nodes are introduced by new couples of compatible labels (e p , e t ) such that e p ∈ L p and e t ∈ L t . Hence, the goal is to sequence as late as possible the labels of L p . To this aim, we first compute the set of labels e t ∈ L t for which the number of labels e p ∈ L p , e p e t is minimal. To break ties, we then choose a label e t such that the average number of labels e t ∈ L t , e p e t , for every label e p ∈ L p , e p e t , is minimal. Then, we introduce in the sequence the selected label e t , preceded by every label e p ∈ L p such that e p e t .

The time complexity of this heuristic algorithm is in

O(n t • log n t • d p • d t ).
Practical computation of an approximate partial order. In practice, one has to compute a total order ≤ i-1 that approximates the partial order i-1 at each iteration i of Algorithm 1. This must be done between lines 7 and 8. Then each call to the test function, line 8, is performed with the total order ≤ i-1 instead of the partial order i-1 . In this case, the time complexity of the computation of i (line 8) is in

O(n t • n p • d • log d).
This complexity can be reduced to O(n t • n p • d) by first sorting all the multisets. When adding the time complexity of the computation of the total order by our heuristic algorithm, we obtain an overall complexity in

O(n t • d • (n p + d t • log n t )).

Filtering within a Branch and Propagate framework

In this section, we introduce two optimizations that may be done when filtering is integrated within a branch and propagate search, where a variable assignment is done at each step of the search.

A first optimization provides an entailment condition for the filtering. If the initial labeling l 0 is such that the maximum label of the pattern graph is smaller or equal to the minimum label of the target graph, every label of nodes of the pattern graph is compatible with all the labels of nodes of the target graph so that no domain can be reduced by our filtering procedure.

A second optimization is done when, during the search, the variable associated with a pattern node is assigned to a target node. In this case, the neighborhood extension procedure is modified by forcing the two nodes to have a same new label which is not compatible with other labels as follows: Definition 11. Let l = (L, , α) be an SIC labeling, and let (u, v) ∈ N p × N t such that v ∈ x u . The propagation of x u = v on l is the new labeling l = (L , , α ) such that -L = L ∪ {l uv } where l uv is a new label such that l uv ∈ L; -= ∪{(l uv , l uv )} so that the new label l uv is not comparable with any other label except itself; α (u) = α (v) = l uv and ∀w ∈ Nodes \ {u, v}, α (w) = α(w)

This labeling l is used as a starting point of a new sequence of labeling extensions. Note that this propagation is done every time a domain is reduced to a singleton.

Experimental Results

Considered instances. We evaluate our approach on graphs that are randomly generated using a power law distribution of degrees P (d = k) = k -λ : this distribution corresponds to scale-free networks which model a wide range of real networks, such as social, Internet, or neural networks [START_REF] Barabasi | Linked: How Everything Is Connected to Everything Else and What It Means[END_REF]. We have made experiments with different values of λ, ranging between 1 and 5, and obtained similar results. Hence, we only report experiments on graphs generated with the standard value λ = 2.5.

We have considered 6 classes of instances, each class containing 20 different instances. For each instance, we first generate a connected target graph which node degrees are bounded between d min and d max . Then, a connected pattern graph is extracted from the target graph by randomly selecting a percentage p n (resp. p e ) of nodes (resp. edges).

All instances of classes A, B, and C are non directed feasible instances that have been generated with d min = 5, d max = 8, and p n = p e = 90%. Target graphs in A (resp. B and C) have 200 (resp. 600 and 1000) nodes.

All instances of class D are directed feasible instances that have been generated with d min = 5, d max = 8, and p n = p e = 90%. Target graphs have 600 nodes. Edges of target graphs have been randomly directed. To solve these directed instances, the filtering procedure is adapted by extending labelings with two multisets that respectively contain labels of successors and predecessors.

All instances of classes E and F are non directed instances that have been generated with d min = 20, d max = 300, and p n = 90%. Target graphs have 300 nodes. Instances of class E are feasible ones that have been generated with p e = 90%. Instances of class F are non feasible ones: for these instances, pattern graphs are extracted from target graphs by randomly selecting 90% of nodes and 90% of edges, but after this extraction, 10% of new edges have been randomly added.

For all experimentations reported below, each run searches for all solutions of an instance.

Comparison of different variants of our filtering algorithm. Algorithm 1 has been implemented in Gecode (http://www.gecode.org), using CP(Graph) and CP(Map) [START_REF] Dooms | Cp(graph): Introducing a graph computation domain in constraint programming[END_REF][START_REF] Deville | Cp(graph+map) for approximate graph matching. 1st International Workshop on Constraint Programming Beyond Finite Integer Domains[END_REF] which provide graph and function domain variables. The global subgraph isomorphism constraint has been combined with c2 constraints (as defined in Section 2.2) and a global AllDiff constraint which are propagated by forward checking.

Table 1 compares different variants of Algorithm 1, obtained by either computing an exact partial order or an approximated one (as described in 4.1 and 4.2), and by considering different limits k on the number of iterations. In all variants, the initial labeling l 0 is the labeling l deg defined in 3.1. Note that the order of l deg is a total order so that in this case the exact and approximated variants are equivalent for k = 1.

Let us first compare the exact and approximated variants. The number of failed nodes with Approx./k = 2 is greater than Exact/k = 2, but it is smaller than with Exact/k = 1. This shows us that the total order computed by our heuristic algorithm is a quite good approximation of the partial order. When considering CPU-times, we note that Approx./k = 2 is significantly quicker than Exact/k = 2.

Table 1 i.e., with the simple l deg labeling: these instances are easy ones, as adding a direction on edges greatly narrows the search space. Instances of classes A, B and C are more difficult ones, as they are not directed; these instances are best solved when k = 1, i.e., after one iteration of the exact labelling extension. Instances of classes E and F, which have significantly higher node degrees, are very difficult ones. For these instances, and more particularly for those of class F which are not feasible ones and which appear to be even more difficult, iterative labeling extensions actually improve the solution process and the best results are obtained when k = 8. As a conclusion, these experimentations show us that (1) the approximated variant offers a good compromise between filtering's strength and time, and (2) the optimal limit k on the number of iterations depends on the difficulty of instances. The best average results are obtained with Approx./k = 4.

Comparison with state-of-the-art approaches. We now compare the variant Approx./k=4 of Algorithm 1 with a state-of-the-art algorithm coming from a C++ library called vflib [START_REF] Cordella | An improved algorithm for matching large graphs[END_REF], and with CP. We consider two different CP models:

-c2 is the model using c2 constraints described in Section 2.2; -c2 + c3 is the model that additionnaly uses c3 constraints introduced in [START_REF] Larrosa | Constraint satisfaction algorithms for graph pattern matching[END_REF].

These two models are combined with a global Alldiff constraint. For c2 and Alldiff constraints, two levels of consistency are considered, i.e., Forward Checking (denoted by FC) and Arc Consistency (denoted by AC). Propagation of c3 follows [START_REF] Larrosa | Constraint satisfaction algorithms for graph pattern matching[END_REF]. All CP models have been implemented in Gecode using CP(Graph) and CP(Map).

Table 2 compares all these approaches and shows us that, except for easy instances of class D which are best solved by vflib, all other classes of instances are best solved by Approx./k=4. When comparing the different CP models, we note that adding redundant c3 constraints significantly improves the solution Table 2. Comparison of state-of-the-art approaches. Each line successively reports the percentage of instances that have been solved within a CPU time limit of 600s on an Intel Xeon 3,06 Ghz with 2Go of RAM; the average run time for the solved instances; and the average number of failed nodes in the search tree for the solved instances.

process except for the easy instances of class D which are better solved with simpler models.

Conclusion

We introduced a new filtering algorithm for the subgraph isomorphism problem that exploits the global structure of the graph in order to achieve a stronger partial consistency. This work extends a filtering algorithm for graph isomorphism [START_REF] Sorlin | A new filtering algorithm for the graph isomorphism problem[END_REF] where a total order defined on some graph property labelling is strengthened until a fixpoint. The extension to subgraph isomorphism has been theorically founded. The order is partial and can also be iterated until a fixpoint. However, using such a partial order is ineffective. Instead, one can map this partial order to a total order. Performing such a mapping is hard, and can be efficiently approximated throught a heuristic algorithm. Experimental results show that our propagators are efficient against state-of-the-art propagators and algorithms. Future work includes the development of dynamic termination criteria for the iterative labeling, the experimental study of other degree distributions, the analysis of alternative initial labelings.

Fig. 1 .

 1 Fig. 1. Instance of subgraph isomorphism problem.

Table 1 .

 1 also shows that the best performing variant differs when considering different classes of instances. Instances of class D are best solved when k = 0, Comparison of different variants of Algorithm 1: Exact (resp. Approx.) refers to the implementation of test described in 4.1 (resp. 4.2); k gives the maximum number of iterations. Each line successively reports the percentage of instances that have been solved within a CPU time limit of 600s on an Intel Xeon 3,06 Ghz with 2Go of RAM; the average run time for the solved instances; and the average number of failed nodes in the search tree for the solved instances.

		Solved instances (%)	Average time	Average failed nodes
		Exact	Approx.	Exact	Approx.	Exact	Approx.
		k=0 k=1 k=2 k=2 k=4 k=8 k=0 k=1 k=2 k=2 k=4 k=8 k=0 k=1 k=2 k=2 k=4 k=8
	A 100 100 100 100 100 100	2.2 0.6 23.4 1.3 1.9 3.2	440 14	0 13	0	0
	B 100 100 100 100 100 100 61.4 5.6 144.2 24.5 28.7 59.6 1314	8	0	3	0	0
	C	45 100 45 100 100 100 439.2 26.3 495.8 101.8 110.4 227.8 1750 13	0	2	0	0
	D 100 100 100 100 100 100 0.7 2.6 99.6 7.5 24.7 56.3	2	0	0	0	0	0
	E	80 60	0 75 80 85 126.7 98.6	-35.2 18.8 36.7 4438 159	-39 13	7
	F	23 20	0 38 63 68 186.4 109.9	-45.0 10.5 3.9 18958 3304	-2323 481 107

If a label e both belongs to Lp and Lt, it is always possible to rename e into e in Lt (where e is a new label), and to add a relation e e for every label e ∈ Lp such that e e.

Acknowledgments

The authors want to thank the anonymous reviewers for the helpful comments. This research is supported by the Walloon Region, project Transmaze (WIST516207).


Laurent Lescop
email: laurent.lescop@nantes.archi.fr

Serge Cassen
email: serge.cassen@univ-nantes.fr

Gavrinis The raising of digital stones

Keywords: Gavrinis, Neolithic, cairn, passage grave, archeology, lasergrammetry, photogrammetry, low-tech techniques, reconstruction, public presentation, real time visits I

The Tumulus of Gavrinis, which is located in Brittany, is one of the best known monuments of Neolithic art because of its extraordinary engravings. A joint team of archaeologists and architects began recording the monument using digital multiscalar technology in 2011. The scanning was complicated by the narrowness of the inner spaces, uncomfortable conditions, granulometry of the stone, and extensive vegetation coverage.

A huge database was created by using a combination of photogrammetric technology and low-tech scanning methods. A lot of post-production work was needed to turn the cloud of data into a clean 3D model, since the original tumulus is irregular and lacks any repetitive patterns or predictive geometry. The goal was to obtain a detailed geometric analysis of the tumulus in order to develop a better understanding of its architecture and engravings, which are sometimes less than a millimeter deep. This article covers the technological solutions that we applied while confronting situations ranging from the need to enhance a carved figure to our effort to reconstruct the entire monument by collating the results of several different recording techniques.

INTRODUCTION

The Tumulus of Gavrinis, on the estuary of the Vannes River in Brittany (Morbihan, France), is known for the extraordinary abundance of engravings on the monolithic pillars that form its walls. The monument has always been a landmark from its erection in the early fourth millennium BC until now. Scientific efforts to describe the engravings have been made ever since the first descriptions in the 19 th century, bringing the site wide scientific and popular renown (Mérimée 1835; Le Roux 1998). Excavations in the seventies and E. Shee-Twohig's precise description of the engravings in 1981 provided researchers with a corpus on which scientific work could be done. The invention of lasergrammetry with millimetric precision made it necessary re-survey the site, since the technology has reached the point where it opens the door to new discoveries.

If the scientific need for a new survey is obvious, the recording must also serve the requirements of specialists in preservation and communication. Gavrinis is highly vulnerable, because of its being one the most famous places in Brittany. The present flow of 30,000 curious and passionate visitors is certainly more than the 16-meter-long by one-meter-wide monument can bear, because of the fragility of its engravings, so better protection must be considered.

There is also a need to celebrate the beauty of Gavrinis' engravings at other sites, in part because a greater appreciation of Neolithic accomplishments can increase the willingness of both the public and decision-makers to protect less famous cairns and standing stones.

Another reason that knowledge of Gavrinis' engravings should be shared with visitors at other contemporaneous sites in the region is that they were originally part of a complex network that has been hidden by rising sea levels, modern construction, and increased vegetation. The landscape that we experience today isolates each Neolithic site as a single entity, but that wasn't the case when they were built.

These communication requirements mean that the new digitized corpus must be suitable for scientific work, use as a training tool for guides, and popularization, both in online applications and in media presentations at a variety of sites. All those different objectives had to be kept in mind during the post-production treatment of the point cloud.

II. OVERALL SCANNING PROCESS

A. Collecting the initial data

The overall process started with a raw Sketchup model based on C.-T. le Roux's and E. Shee-Twohig's renderings. As Gavrinis is situated on an island, which is only accessible by boat, and has a poor electricity supply, everything had to be carefully prepared. The main problem that emerged was the narrowness of the interior spaces. It forced us to imagine the scanning process as a progression from outside to inside, with a large overlap.

The survey involved recording geolocational data with different scales of acquisition. Thus the whole cairn and surrounding ground were scanned using a Leica Geosystems C10 which is a time-of-flight system. Then, the walls of the tomb and the upper face of the cover slab were recorded with a Nikon Krypton K610 scanner with infra-millimetric resolution associated with Kube software.

For the Leica recording, six stations were set around the cairn, white spheres were used as targets to help the merging process. Eight more stations turned out to be necessary inside the corridor to link the inside to the outside. This system was selected for its polyvalence, allowing a single campaign to record both the architecture and slabs. The Leica system has a range of 300 m (albedo 0.9), 134 m (albedo 0.18) and a minimum range of 0.1 m. The resolution depends to the laser beam which, in this case, was between 0 to 50 m: 4.5 mm (based on FWHH), 7 mm (based on Gauss).

The orthostats in the passage are registered with a letter L for left and R for right and a letter based on Shee-Twohig's nomenclature. They were recorded with a Krypton K610 from Nikon Metrology at an infra millimetric scale. This scanner works with a 200-mm-wide laser beam. Its potential precision is 37 µm but a precision level of 0.5 mm was considered sufficient, for the moment. It turned out that this level wasn't precise enough for further analysis, but the gap was filled by using by an analogic oblique lighting method [START_REF] Cassen | forthcoming. In search of a representation of the Gavrinis megalithic tomb (Morbihan, France), à paraître[END_REF]

and forthcoming).

Each engraved pillar holds three to five million points. Focusing software was used to clean up the data and transcribe it into both XYZ and STL formats. The recording of the pillars represents 72,399,271 points. The recording of the cairn and its surroundings required 3,937,900 points. After eliminating redundant information and extraneous features such as plants and trees, we got a total of 93,564,668 points geo-recorded in the Lambert II and IGN 69 coordinate systems. The merging of all the point clouds was done in Leica Cyclone. The IMP files were then exported in XYZ, PTS, and TXT. A first raw meshing was exported in STL and OBJ.

B. First results

The first result provided a precise cross-section of the cairn, positioning the corridor very precisely in the mass of the monument. From the point cloud, it was quick and easy to draw a ground plan and elevations as well as horizontal and vertical sections. Unlike many other passage graves (Mané Lud, Mané Rutual, Petit Mont, Table des Marchands), no precise blueprints existed for Gavrinis. An angular deviance of 25 cm exists between the last rendering (Le Roux 1985) and the one we have now in the path from the the entrance to the back of the chamber.

The cross-section of the passage confirmed the presence of a previous mound underneath the one we know (Closmadeuc 1864; Cassen, [START_REF] Boujot | La Constitution sédimentaire des tertres[END_REF]. The angle of the slope, can be measured (2°) even if we know the ground has been excavated and may be a bit different. The laser recordings also made it possible to measure and compare the depths of the engravings. Some seemed to be more altered in a certain direction, confirming the hypothesis of former use before the engraved blocks were in the cairn (Le Roux 1984). The erosion, visible with a simple flashlight, cannot have happened inside the enclosed space; the stones must have been standing outside before being embedded in the cairn. The atmospheric effects, meteorisation, are now highly visible and can be shown in sectional views.

III. ENHANCING THE ENGRAVINGS

Twenty-three out of the 26 pillars are engraved with a complex network of lines, figurative objects and dense patterns. Motifs cover the entire surface of pillars, often slipping around the sides. The details show the astonishing richness of the composition, which plays with the initial shape of the stone.

A. The L6 case

We will use Pillar n°L6 (Meaning the one in the sixth position from the entrance on the left) as our first example because of its iconographic and recording challenges and enormous potential. The purpose here is not to describe the precise significance of signs and patterns; but to depict the means that led to better analyses.

The first idea was to mimic an analogical process with software. Using oblique lighting is highly efficient in revealing details, but it takes time, requires being on location and is sometimes uncomfortable, because of the size and shape of the the spaces being recorded. A CG process on a dense geometry would give the same results.

The first tests were made with Meshlab, but despite of the free software's benefits, the processing took too long: every step took dozens of minutes to complete. VrMesh was also tested, but the processing was still long. This software now specializes in Lidar data recognition. In the end, Geomagic was chosen after several tests, both because it finished its calculations faster than the competing programs and because it offered a lot of tools for building perfect geometries, which was our goal for further real-time use.

A mesh of 2,501,332 triangles was formed out of the point cloud and virtual lighting was used to simulate oblique lighting. The flexibility offered by this approach proved to have added value over the analogous method. A further refinement was added to complement these results. The geometry was copied and then smoothed to a high degree. This new geometry is intercalated with the original one, in a new color, filling in the lowest parts of the engravings and giving the whole a flatter look, which allows the engravings to be more easily read in relationship with each other. The idea of having an intersection layer was presented by J.A. Sanz [START_REF] Sanz | A simple methodology for recording petroglyphs using low-cost digital image correlation photogrammetry and consumer-grade digital cameras[END_REF], who demonstrated its utility with a single sample, whose surface was flatter than ours. A relatively flat shape was placed virtually halfway down the engravings, isolating their upper parts. Our technique is very similar, but has to work with more complex 3D surfaces. It requires step by step adjustments and a lot of fine tuning, depending to the depths of the engravings, granulometry of the stone, and overall shape of the orthostats. Many tests were necessary to find the optimal smoothness for achieving the greatest clarity. Both the original and smoothed geometries are shown in our results.

The deviation between the smoothed and original surfaces also produced good results, which allowed us to extract further information and provide a better synoptic view of the engravings. For example, it was possible to use the color spectrum to make the depth and width of all the lines within a 6 mm range clearer than they were in the original renderings. This work was done in Geomagic. Interesting details then emerged in this deviance analysis. Some figures became more recognizable and could be found carved on several standing stones. One of them, with a triangular shape, is an axe, replicating the exact form of a polished axe. The turning light technique (which consists of taking a large series of photos with oblique light oriented at a different angle for each view), detected a very interesting detail that could be isolated and measured. A tiny hole in the figure of an axe, previously barely visible, looked, in size and shape, to be very similar to those observable on real axes [START_REF] Pétrequin | Chronology and distribution of large Alpine axeheads in western Europe, in Typology, chronology and distribution of large Alpine axeheads, Presses Universitaires de Franche-Comté n°1224, Collection Les cahiers de la MSHE Ledoux n°17[END_REF] When the zone with the triangular shapes was mapped topographically, the curves formed by the contour lines revealed the holes. This procedure provided a further way of detecting and measuring an interesting element. Further analyses of these details are the domain of specialists, among them, traceologists. The first step in interpreting the engravings is establishing their chronology. Knowing their temporal sequence provides clues for understanding what we will provisionally call their "nonfigurative signs". Making chronological assessments requires the use of both real oblique lighting and the digital techniques described above. The overall process can be summarized in the following drawing: And it has to be done for each pillar.

B. The covering slab

The slab that covers the chamber is well known for being part of a gigantic stone broken into 3 main pieces, which were reused in two major monuments. The figurative engravings on the upper side of the slab fit ones at La Table des Marchands. The theory that the slabs at the two sites were originally joined is confirmed by the fact that their broken edges match each other perfectly.

The covering slab is only visible in the chamber and from above, in a space that was made over the chamber in modern times for viewing the engravings. The laser scanning captured both the top and bottom sides, but, of course, the embedded parts were inaccessible to the scanner. The stone itself is therefore only visible in the resulting renderings only as a vacuum. By isolating the two surfaces, it becomes possible to understand how the slab lies on the pillar. The thickness is also measurable. The two sides were imported into Mudbox for overall reconstruction. Starting from a basic volume, the sides were applied like stencils, transferring all their tiny details to the new geometry. The invisible parts were then sculpted. The upper part of the covering slab is only visible from a room created during the monument's restoration in the early seventies. The fact that the floor of this technical space is precisely the size of the engraved slab makes it impossible to put the laser scanner anywhere else than on the slab.

The narrowness of the chimney, which provides access to the space via a ladder, makes it impossible to use many scanning systems. It even made it difficult to get the scanner, cameras and other small devices into the space and impaired our ability to completely record the slab. As a consequence, large parts are absent from our laser scan. Hopefully, they are not among the most interesting areas.

We decided to use photomodelling to complete the model and add color information that didn't exist in the point cloud. A Canon EOS 600D camera with a Canon EF-S 18-55mm f/3.5-5.6 IS II lens was used for this purpose. Two projectors with natural daylight colors were positioned to increase the light and distribute it in the most homogeneous way.

The main problem in such an enclosed and dark place is to avoid projecting one's own shadow upon the subject of the photography, especially when it is on the floor. The overlapping of views also has to be anticipated. In this case, it was necessary for the photographer to revolve in place over a specific detail.

Photoscan was used. The model was completed using 500,000 vertices, which produced 250,194 polygons. The texture was calculated as an orthophoto which can be used as a texture. The laser model and orthophoto obtained by photomodeling were blended to fill in areas missed by the laser and provide texture. A clean texture with proper UVW coordinates can be used for a real-time model on which it is possible to append graphic information to make it easier to interpret. This side of the slab is not open to visits by anybody except researchers. The only representations are drawings in books. Our results make it possible to examine the engravings on the upper side of the slab for the first time by visiting the space above them virtually. C. The threshold slab.

Another step in making an overall virtual model of Gavrinis was scanning parts of the tumulus, which are now held elsewhere, so they could be re-incorporated in the final digital model.

The threshold slab, which is registered as S12, is kept in a repository of the Regional Archaeological Service of Brittany in Vannes. Three recordings were carried out on this specific stone. The first used the turning light method, which emphasizes details; the second was a regular scanning process with the same device, the Krypton K610, used in the cairn; the third one was photomodeling.

The resolution chosen for the scanner was 0.5 mm. After the cleaning and meshing process, we got a total of 17,274,685 polygons. The deviance technique was used to enhance the details and better define the engravings. The third recording for photomodeling involved photogrammetry. We used a Canon EOS 500D at 18 mm. The subject was exposed in an ideal diffuse light environment. 14 photos were taken for the top side and then 29, for the rest of it. Visual SFM and 123Dcatch were then compared for the reconstruction and both of the results to the laser recording. 123D Catch is a free online application offered by Autodesk. It is a "black box", meaning that it only has a few tweaking solutions, but the reconstruction is usually good. In this case, the deviation of 123DCatch was only 1.78mm.

Visual SFM is software that aggregates the photo-matching process, point cloud generation, and meshing. A lot of finetuning possibilities make it possible to produce a good survey. VSfm only deviated in this case by 0.98mm. The deviance between the laser scan and photomodeling was not uniform. Depending to the light and number of photos used, the differences were between one and two centimeters, which provided satisfactory results for this kind of subject. 123Dcatch, which is free software, is advertised as offering real scale models. A short test showed that the results were close to the reliable ones provided by lasergrammetry. The problem remains that cameras are hardly able to reveal geometric details except if photos are taken at a very short distance for a limited surface or, if the light is oblique enough. Tests therefore should be made to measure the depth of hollows and their orientation.

D. P2 The hidden engraving.

Recording surfaces in tight spaces cannot be done without proper knowledge of the place. It requires an overall strategy, involving decisions as to where to put devices, what kind of moves to make, what parts need to be focused on, and what results one is aiming for.

We knew for sure that an important sign was carved in Gavrinis's corridor. Due to close proximity between the motif and next slab, it was very difficult to record it, and we could not find a device that could slip into the gap. The widest part is only 14cm, and then it quickly shrinks to 1cm in a 40cm deep interstice. The cairn is a dome, with a diameter of 25m and height of 8m. The main façade, which has been reconstructed, is a 30m long stone wall with grading. The rest of the monument is covered by vegetation, most of which is weeds. Before any reconstruction, the first task was to "mow" the tumulus by virtually eliminating all the plants and other unwanted natural information. Holes in the point cloud were filled. The stone facade offers lots of details of different size and automatic construction processes can be thwarted or generate highly degraded geometries showing spikes, tunnels, selfintersections, creased edges or small holes, even though the eye easily catches the overall shape and structure. The meshing was re-arranged and simplified (the goal with to stick between 50,000 to 300,000 polygones) to increase the regularity, aid volume calculations and make virtual visits in real time easier. This adjusted topographical model was in keeping with the details of a normal map and a depth map can easily be added.

Geomagic offers tools for reaching cleaner geometries, which are essential for any further post-production treatments. The retopology stage can be performed by many software programs, but it must often be carried out by hand and can be very time-consuming. The latest version of Mudbox provides another nice and easy tool for adjusting the topological parameters and the learning curve is acceptable. Using the model for real-time applications becomes easy. Gavrinis is a complex subject and more than one tool was needed to achieve our goals, which were to record, interpret, speculate and communicate. The overall digital reconstruction of Gavrinis is still in progress. Our work raised new problems, techniques and questions as it progressed. Neolithic funeral architecture is intractable: hidden parts often cannot be deduced from what can be seen, each part is slightly different from all the others, and yet the monuments are laid out precisely. Whether lasergammetry or photogrammetry is used, the recording process is usually constrained to confined spaces, obscurity and long distances from power resources.

That's the reason why we added a variety of lighting and low-tech methods to our toolkit. Over our years of practice, which has teamed archeologists with architects, we have become deeply convinced of the importance of using several complementary devices and tools on the same subject. We are not just digging for information or trying to produce a better 3D computer model, but are trying to articulate scientific needs with the requirements of satisfying visitors' curiosity.

Fig. 1 .

 1 Fig. 1. Gavrinis, outside view of the cairn.

Fig. 2 .

 2 Fig. 2. Krypton K610 from Nikon Metrology associated with Kube software.

Fig. 3 .

 3 Fig. 3. The position of the corridor in relation to the mass of the cairn.

Fig. 4 .

 4 Fig. 4. The assembled point cloud of Gavrinis.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Groundplan of Gavrinis' corridor

Fig. 7 .

 7 Fig. 7. Exemples of orthostats : L7 and L5

Fig. 8 .

 8 Fig. 8. L6, digital oblique lighting tests.

Fig. 9 .

 9 Fig. 9. L6, ungravings revealed by a smooth intercalated geometry.

Fig. 10 .

 10 Fig. 10. L6, deviance analysis from the original mesh and a smoothed one.

Fig. 11 .

 11 Fig. 11. L6, deviance showing hole details in the figure of an axe. A third hole can be seen.

Fig. 12 .

 12 Fig. 12. L6, hole corresponding to real objects.

Fig. 13 .

 13 Fig. 13. L6, analysis process for one stone.

Fig. 14 .

 14 Fig. 14. The connection of slabs as a single piece (Cassen).

Fig. 15 .

 15 Fig. 15. The slab is revealed as a vacuum between two surfaces.

Fig. 16 .

 16 Fig. 16. P10, the two surfaces and view from above.

Fig. 17 .

 17 Fig. 17. P10, The yellow parts shows the theoretical reconstruction of inaccessible embedded portions of the slab.

Fig. 18 .

 18 Fig. 18. P10, orthophoto of the upper face.

Fig. 19 .

 19 Fig. 19. P10, CG model for real-time use with interpretation.

Fig. 20 .

 20 Fig. 20. S12, Mesh model from K610 laser scanner .

Fig. 21 .

 21 Fig. 21. S12, samples of photos used for photomodeling.

Fig. 22 .

 22 Fig. 22. Deviation between laser and 123Dcatch .

Fig. 23 .

 23 Fig. 23. Deviation between laser and VSFM.

Fig. 24 .

 24 Fig. 24. P2, position of the engraving We finally used a Nikon D5000 at 18mm to get a collection of 32 photos (RAW format 4288x2848) and PhotoScan for the reconstruction. After cleaning, we had 72,471 polygones. It was possible to get a proper elevation and a front projection.

Fig. 25 .

 25 Fig. 25. P2, reconstruction with PhotoScan and front projection

Fig. 26 .

 26 Fig. 26. Cairn, huge mixture of vegetation and architecture

Fig. 27 .

 27 Fig. 27. Cairn, before and after retopology, then textured for real time use.

TABLE I .

 I

	LASERGRAMMETRY AND PHOTOMODELING COMPARAISON
	Technique	Length	Width	Thickness	Volume	Deviance
		mm	mm	mm	m 3	mm
	Lasergrammetry	1274	794	296	0.218	-
	123DCatch	1272	797	298	0.225	1.78
	VSFM	1279	814	296	-	0.98

TABLE II .

 II LASERGRAMMETRY AND PHOTOMODELING COMPARAISON

	Technique	Length	Width	Thickness	Volume	Mass
		mm	mm	mm	m 3	kg
	Lasergrammetry	1274	794	296	0.218	589
	123DCatch	1272	797	298	0.225	607

ACKNOWLEDGMENT

The authors wish to thank the Ministère de la Culture et de la Communication, le CIRA (Commission interrégionale de la recherche archéologique) le CG56 (Conseil Général), le SDA (Sous Direction de l'Archéologie). We are also very grateful to Y.Bellenfant, SAGEMOR (Société Anonyme de Gestion du Morbihan), the CNRS (Centre national de la recherche scientifique) and ENSANantes (Ecole Nationale Supérieure d'Architecture de Nantes). The authors also wish to warmly thank Duncan Caldwell for revising the text.


Antoine Muller
email: antoine.muller@irisa.fr

Félix Demore

Charles Pontonnier

Georges Dumont

MusIC MAKES THE MUSCLES WORK TOGETHER

INTRODUCTION

In most musculoskeletal simulation software, the muscle force estimation step is solved thanks to an optimization problem. It consists in minimizing a criterion representing the central nervous system strategy. Despite several improvements, this optimization remains computationally costly. Some authors propose real-time simulations but achieve compromises on model complexity [3] or on the optimization stop criterion [1] to achieve this performance.

In this work, we present a fast and quasioptimal method of muscle forces estimation: the MusIC method (Muscle forces Interpolation and Correction). It is based on precomputed data collected in a database as previously proposed by [START_REF] Pontonnier | Symposium on CMBBE[END_REF]. Its performance in terms of computation frequency and accuracy is showed on a usually used musculoskeletal leg model thanks to a cycling motion.

METHODS

From joint coordinates q and joint torques Γ, the MusIC method aims at computing a muscle forces vector solution F. It is composed on two steps.

The first one consists in generating a database beforehand by using only the musculoskeletal model (Database generation in Fig. 1). Considering each joint separately, it contains an activation ratio vector, solution of the classical optimization problem for several joint configurations, ranging the physiologically feasible joint space. This vector corresponds to the activation part of each muscle compared to the sum of activations of all these muscles.

In the second step, from the current state, we use the database to find a solution of the force sharing problem for each joint (Interpolation in Fig. 1). The activation ratio vector is extracted from the database thanks to a linear interpolation with the current joint configuration. The solutions of each joint are then corrected to compute a muscle forces vector satisfying the global dynamic equilibrium and the physiological properties (Correction in Fig. 1). For this, a single muscle forces vector is computed thanks to a barycentric interpolation with regard to the different joints. The closest solution of this muscle forces vector which satisfies the global dynamic equilibrium and the physiological properties is then considered as the muscle forces vector solution.

This method was compared to a classical optimization problem solved thanks to a SQP method. For each of them, a 3-order polynomial criterion is applied where each muscle force is normalized by its maximal isometric force. One male participant (184 cm, 70 kg) performed a cycling motion. 47 motion capture markers were placed on standardized anatomical landmarks and captured thanks to a Vicon ® motion capture system (125 Hz). Two instrumented pedals (1000 Hz) were used to access external forces applied on each foot. We reconstruct the motion thanks to a whole body model. Musculoskeletal model of the legs is extracted from the work of [2] and is composed of 82 muscles distributed over the two legs.

RESULTS AND DISCUSSION

From muscle forces solution of both methods, we recompute the cost function value throughout the motion. The relative difference between the average values is 2.7%. Thus, the solution suggested by the MusIC method is close to the classical one in terms of optimality. Finally, the computation frequencies are respectively 1.7 Hz and 36 Hz for the classical optimization problem and for the MusIC method. As the database generation is performed beforehand, the computation time of this step is not taken into account.

The results described above show that the MusIC method finds a quasi-optimal solution. Its main benefit is the computation frequency which is approximately twenty times higher. One of the principles of this method is to consider separately the muscles contribution for each joint. However, despite the significant part of poly-articular muscles, the solution seems consistent with the optimal one. Moreover, the use of a linear interpolation seems to filter the muscle forces compared to those obtained with a classical optimization. For example, we can see in Fig. 2 that the rectus femoris force curve obtained with the classical method exhibits an artefact. It comes from a local minimum in the optimization method. As the MusIC method is based on joint configurations which were filtered throughout the motion, this artefact does not appear with this method.

The MusIC method is based on a pre-computed database which depends only on the musculoskeletal model. Thus, it is usable whatever the motion of a subject is and can improve largely the computation time to achieve real time computation. This method could be extended to other cost functions, for example considering the joint reaction forces. Obviously, additional validations have to be performed with other musculoskeletal models to completely validate the MusIC method.

CONCLUSIONS

We presented a fast and quasi-optimal method of muscle forces estimation: the MusIC method, which was partially validated on a musculoskeletal leg model. The solution obtained is close to the optimal one with a particularly low computation time compared to a classical optimization method.

Fig 2 :

 2 Fig 2: Muscles forces (each normalized by its maximum isometric force) on the right leg throughout the motion. The continuous blue line and the dotted red line correspond respectively to results obtained with the classical optimization problem and with the MusIC method.

Fig. 2

 2 Fig. 2 shows representative examples of muscle forces (each normalized by its maximum isometric force) obtained with both methods. For each of them, periods of activation and non-activation and amplitude of activation are consistent. Considering all muscles, the mean correlation coefficient of the MusIC method results is 0.89(±0.14) compared to the classical optimization problem.