Hélène Touzet
email: helene.touzet@univ-lille1.fr

On the Levenshtein Automaton and the Size of the Neighbourhood of a Word

Given a word P and a maximal number of errors k, we address the problem of counting the number of strings whose Levenshtein distance to P does not exceed k. We give an algorithm that scales linearly with the size of P and that is based on a variant of the classical Levenshtein automaton.

Definition and Construction of the Deterministic

Universal Levenshtein Automaton

Preliminaries and notations

Let Σ be a finite alphabet, and let P and V be two strings over Σ. The Levenshtein distance between P and V , denoted Lev(P, V ), is the smallest number of edit

Introduction

The problem of measuring the similarity between two strings arises in many areas such as computational molecular biology, natural language processing, spelling correction, plagiarism detection, music information retrieval. A common metric for it is the Levenshtein distance, also simply called the Edit Distance. This distance is defined as the smallest number of substitutions, insertions, and deletions of symbols required to transform one of the words into the other.

In this paper, we investigate the basic problem of the size of the neighbourhood of a given string: count how many strings are within a bounded distance of a fixed reference string. This problem has been exposed in [START_REF] Myers | What's behind blast[END_REF][START_REF] Becerra-Bonache | Learning balls of strings from edit corrections[END_REF], among others. As far as we know, there is no efficient algorithm for solving it. We propose a dynamic programming algorithm that runs linearly in the length of the pattern. This algorithm heavily relies upon the Universal Levenshtein Automaton, which is an advanced automaton for the Levenshtein distance problem introduced in [START_REF] Mihov | Fast approximate search in large dictionaries[END_REF][START_REF] Mitankin | Universal Levenshtein automata. Building and properties[END_REF]. This automaton is universal in the sense that it does not depend on the two strings that have to be compared.

The paper is organized as follows. In Section 2, we present the Universal Levenshtein Automaton. We review its main principles, and revisit them by proposing some new features, such as the introduction of a Nondeterministic Universal Levenshtein Automaton. In Section 3, we present the algorithm to compute the cardinality of the neighbourhood of a word. operations needed to transform P into V , where the allowed operations are substitution of one symbol with another, deletion of a symbol and insertion of a symbol. A sequence of operations transforming P into V is called an edit script. For any pair of strings P and V , the distance between P and V can be computed by dynamic programming in time proportional to the product of the length of P and V with the Wagner-Fischer algorithm [START_REF] Wagner | The string-to-string correction problem[END_REF].

Assume now that we are given a fixed threshold for the number of errors k. We consider the decision problem associated to the Levenshtein distance: for any pair of strings P and V , decide whether the Levenshtein distance is lesser than or equal to k. In this decision problem, we can also introduce the additional hypothesis that the pattern P is fixed, and consider the neighbourhood of P , noted Lev(P, k):

Lev(P, k) = {V ∈ Σ * ; Lev(P, V ) ≤ k}.
When P is fixed, it is highly desirable to preprocess P , so that the decision problem can be solved more efficiently than with the dynamic programming algorithm. There has been an abundance of literature on this subject. A standard approach is to start from the Levenshtein automaton, depicted in Figure 1, which recognizes the set of strings which are at most at distance k to P . This Levenshtein automaton is a Nondeterministic Finite Automaton (NFA). A run can be simulated using dynamic programming or bit parallelism [START_REF] Baeza-Yates | A faster algorithm for approximate string matching[END_REF][START_REF] Holub | Implementation of nondeterministic finite automata for approximate pattern matching[END_REF]. However, it is of little help to count the number of accepted strings, because distinct paths in this NFA can recognize the same string. Another possibility is to transform the Levenshtein automaton into an equivalent Deterministic Finite Automaton (DFA) [START_REF] Ukkonen | Finding approximate patterns in strings[END_REF][START_REF] Schulz | Fast string correction with levenshtein automata[END_REF]. This is a tedious task, that should be performed for each target string P .

In [START_REF] Mihov | Fast approximate search in large dictionaries[END_REF][START_REF] Mitankin | Universal Levenshtein automata. Building and properties[END_REF], a universal Levenshtein automaton was introduced. The term universal conveys its one-time construction and independency of the two input strings. Thus it can be applied to any pair of strings P and V of arbitrary length over any arbitrary alphabet Σ. This new automaton is based upon insightful observations of the nondeterministic Levenshtein automaton for a fixed word, and extends the work done for the determinization algorithm described by the same authors in [START_REF] Schulz | Fast string correction with levenshtein automata[END_REF].

The remainder of this section is devoted to a thorough presentation of the Universal Levenshtein Automaton. First, in Section 2.2, we explain how to convert the two input strings, P and V , into a single string that contains the required information to determine their distance. Then, in Subsections 2.3 and 2.4, we explain how to construct the Universal Levenshtein Automaton.

Notation: For a word V ∈ Σ * , |V | is the length of this word, and

|V | Σ is the number of distinct symbols of Σ present in V . For each positions i, j in V (1 ≤ i ≤ |V |), V i is the ith letter of V and V [i.
.j] is the subtring of V starting at position i and ending at position j. ε denotes the empty string.

Bit vector representation

Let P and V be the two strings. In the two-dimensional dynamic programming table of the Wagner-Fischer algorithm, all edit scripts with at most k errors have

S S S A A A L L L A A D D D A Σ Σ Σ Σ Σ Σ Σ Σ Σ Σ Σ Σ Σ , ε Σ , ε Σ , ε Σ , ε Σ , ε Σ , ε Σ , ε Σ , ε Σ , ε Σ , ε 0 #0 1 #0 2 #0 3 #0 0 #1 0 #2 4 #0 1 #2 2 #2 3 #2 4 #2 5 #2 1 #1 2 #1 3 #1 4 #1 5 #1 5 #0
Fig. 1. Levenshtein automaton for the word SALAD. In general, this automaton has m + 1 columns and k + 1 levels. The i #x notation for each state corresponds to i symbols read in the pattern and x errors recorded. Horizontal transitions represent identities, vertical transitions represent insertions, and the two types of diagonal transitions represent substitutions (Σ) and deletions (ε), respectively.

a path that stays around the diagonal. The width of this diagonal is 2k + 1. To compute the values in this portion of the table, it is sufficient to know which coordinates (i, j) of the cells are such that P i = V j . This information captures the similarity between the two strings P and V . We show an example in Figure 2 (left), where grey cells represent coordinates such that P i = V j , and white cells represent other coordinates. At this point, one can forget the two input words, P and V . The idea of the Universal Levenshtein Automaton is to work directly on these configurations of white and grey cells. For this, each horizontal line of the diagonal is represented by a bit vector of length 2k + 1: 0 for white cells, and 1 for grey cells. These bit vectors will serve as the new alphabet, rather than Σ. This encoding captures the local structural properties of the input words, and guarantees alphabet independence.

Definition 2. Let P ∈ Σ m , V ∈ Σ n such that n ≤ m + k. Let P = $ k P $ 2k and V = V $ m-n+k
, where $ is a new symbol not present in Σ. The k-encoding of V with respect to P is the sequence of m + k bit vectors of length 2k + 1 such that the jth element is the characteristic vector χ(V j , P [j -k..j + k]).

The new symbol $ is a sentinel character that serves two purposes. It is added to the prefix and to the suffix of P to standardize the length of bit vectors. Additionnaly, it is added to the suffix of V to deal with edit scripts that end with one or several deletions of symbols of P . These operations will be treated as substitutions with $. Note that the k-encoding of V with respect to P is not defined when the length of V exceeds the length of P by more than k characters. Indeed, there is no point in asking Lev(P, V ) ≤ k in this case. Finally, we define Lev(P, k) as the set of bit vector sequences u in ({0, 1} 2k+1 ) * such that there exists V in Lev(P, k) whose k-encoding wrt P is u. Remark 4. In [START_REF] Mihov | Fast approximate search in large dictionaries[END_REF], the authors use a different encoding. They have bit vectors of length 2k + 2 bits, instead of 2k + 1. The last bit is used to identify the transition between non-accepting and accepting states. Here, we get rid of this additional bit by adding $ symbols to the suffix of V . In Subsection 2.4, it will allow us to obtain a smaller automaton.

Construction of the Nondeterministic Universal Levenshtein Automaton

In [START_REF] Mihov | Fast approximate search in large dictionaries[END_REF], the authors directly build the Deterministic Universal Levenshtein Automaton from the Levenshtein automaton for a fixed word. They consider a symbolic triangular area of a state to simulate multiple active states. Here we take a different approach, and introduce a Nondeterministic Universal Levenshtein Automaton. We find this intermediate step useful to facilitate understanding. Moreover it allows us to give an effective algorithm to build the Deterministic Universal Levenshtein Automaton, which is not so easy to infer from [4] 1 .

Let us come back to the table in Figure 2. An edit script between P and V can be seen as a path in the white and grey grid. The portions of the path that stay in the same lane, correspond to a series of identity and substitution operations. In this context, traversing a white cell costs one error. Everytime you change lanes, you have to pay either for an insertion (moving to the left) or a 1 The authors write: "We describe only the basic idea." deletion (moving to the right). To capture this idea, we introduce states of the form (x, y), meaning "I am in the lane y, and have made x errors so far".

We now examine which are the outgoing transitions for the state (x, y). If we prefer not to have ε-transitions, the automaton should read a bit vector from its input sequence at each time step. So each transition should consume exactly one symbol v of V . In a first approach, an edit script can be decomposed into a series of one the two following basic events: an insertion of v, or a series of deletions in P (0 ≤ ≤ k -x), followed by either a substitution or an identity with v. In this decomposition, we make the usual hypothesis that a deletion is not followed by an insertion. To reduce the nondeterminism, we add another local condition: a deletion is not followed by a substitution. Indeed, it is always possible to intervert the two operations, and to apply the substitution before the deletion. Doing so, we obtain three types of edit events.

ins: insertion of v, -sub: substitution of v, -del+id: deletions in P (0 ≤ ≤ k), followed by an identity with v. ins makes the automaton transits from (x, y) to (x + 1, y -1), sub from (x, y) to (x + 1, y) and del+id from (x, y) to (x + , y + ). Figure 3 shows all these transitions.

-0 1- -0 - -0 - -0 01- 2 del+id 3 del+id 0 del+id -0 001- -1 - 1 del+id sub x, y x + 1, y x + 2, y -1 x + 3, y + 3 y + 2 x + 1, ins y + 1
x + 1, Fig. 3. Outgoing transitions for the state (x, y), k = 3. For each bit vector, the bit in position k + y + 1 corresponding to the lane y is framed, and -is used to denote any sequence of bits.

We are now ready to formally define NULA(k), the Nondeterministic Universal Levenshtein Automaton for k errors. the input alphabet is {0, 1} 2k+1 , the set of states

Q k is {(x, y) ∈ IN × Z Z; 0 ≤ x ≤ k, -x ≤ y ≤ x}, -the transition function ∆ k : Q k × {0, 1} 2k+1 → P (Q k ) is constituted of three types of transitions.
insertion transitions: (x + 1, y -1) ∈ ∆ k ((x, y), u) for all states (x, y) ∈ Q k such that x < k and all u ∈ {0, 1} k+y 0{0, 1} k-y substitution transitions: (x + 1, y) ∈ ∆ k ((x, y), u) for all states (x, y) ∈ Q k , such that x < k and all u ∈ {0, 1} k+y 0{0, 1} k-y deletion+identity transitions: (x + , y + ) ∈ ∆ k ((x, y), u) for all states (x, y) ∈ Q k , all such that 0 ≤ ≤ k -x and all u ∈ {0, 1} k+y 0 1{0, 1} k-y- the start state is (0, 0), all states are accepting.

Figure 4 shows NULA(2). The automaton is very regular and easy to construct. The reader can check that all 2-encodings of Example 3 are accepted. In general, NULA(k) can be seen as an extension of NULA(k -1): just add 2k + 1 new states of the form (k, y) and incoming transitions. Thus NULA(k) has (k + 1) 2 states.

•1••• 1•••• ••••1 ••001 ••0•• ••01• •••0• ••0•• ••01• ••1•• •1••• ••1•• •••1• •••1• •01•• •0••• •0••• ••0•• ••0•• •••0• 1,0 0,0 2,-2 2,- 1 2,0 2,2 1,-1 2,1
•••01 1,1 ••1•• Fig. 4. NULA(2)
, the Nondeterministic Universal Levenshtein Automaton for k = 2 errors. All states on the same horizontal level carry the same number of errors x, and all states in same same column correspond to the same phase y in the pattern. Each looping arrow is an identity transition, each vertical arrow a substitution transition, each north-east arrow a deletion transition, and each north-west arrow an insertion transition. For the sake of readability, the symbol • in a bit vector is either 0 or 1. (0, 0) is the start site and all states are accepting.

We prove that NULA(k) effectively recognizes the expected language. Recall that the right language of a state q ∈ Q k , denoted L(q), is the set of all sequences of bit vectors u such that NULA(k) when started in q will accept u.

Proposition 6. Let (x, y) ∈ Q k . Given P ∈ Σ m and V ∈ Σ n , we have when y = 0: u ∈ L(x, y) if, and only if, Lev(P, V ) ≤ k -x, when y > 0: u ∈ L(x, y) if, and only if, Lev(P [1 + y..m], V ) ≤ k -x, when y < 0: u ∈ L(x, y) if, and only if, Lev(P,

V [1 -y..n]) ≤ k -x,
where u is the k-encoding of V with respect to P .

Proof. The proof is by induction on the length of u.

Corollary 7. Let P ∈ Σ * , V ∈ Σ * . Lev(P, V ) ≤ k, if, and only if, the kencoding of V with respect to P is accepted by NULA(k).

Construction of the Deterministic Universal Levenshtein Automaton

We now explain how to build the Deterministic Universal Levenshtein Automaton from the NFA introduced in the preceding section. The principle is to use the standard powerset construction that converts a NFA into a DFA. In this construction, we show that it is possible to reduce the number of states by defining subsumed states. This notion is already present in [START_REF] Schulz | Fast string correction with levenshtein automata[END_REF], where it is defined for the Levenshtein automaton for a fixed word. Here, we adapt it to the states of NULA(k).

Definition 8. Let (x, y) and (x , y ) be two states of Q k . We say that (x, y) subsumes (x , y ), denoted (x , y ) (x, y), if x < x and y+x-x ≤ y ≤ y+x -x.

It is clear from the definition that the relation is a well-founded partial order.

Proposition 9. Let q ∈ Q k , q ∈ Q k , such that q q. Then L(q ) ⊆ L(q).

Proof. Consequence of Proposition 6.

This proposition implies that all subsumed states can be removed from a subset of Q k without modifying its right language. In other words, for any subset Q of Q k , Q and Reduced(Q ) cannot be distinguished, where Reduced(Q ) is defined as the largest subset of Q such that no two elements of Reduced(Q ) are subsumed. It allows us to prune the set of states considered during the construction of DULA(k). Figure 1 gives the corresponding algorithm and Definition 10 the formal definition of DULA(k).

Definition 10. Let k be a positive number. The Deterministic Universal Levenshtein Automaton for k, denoted DULA(k), is the DFA represented as follows.

the input alphabet is {0, 1} 2k+1 , the set of states is the set of reduced subsets of Q k , the transition function δ is given by Algorithm of Figure 1, the start site is {(0, 0)}, all states are accepting.

Add {(0, 0)} to DU LA(k) as an unmarked state; while DU LA(k) contains an unmarked state do Let T be that unmarked state; Mark T ; for each bit vector u ∈ {0, 1} 2k+1 do S = {q ∈ Q k ; ∃q ∈ T q ∈ ∆(q, u)}; S = Reduced(S); Define δ(T, u) = S ; if S is not in DU LA(k) already then Add S to DU LA(k) as an unmarked state; end end end Algorithm 1: Construction of DULA(k) from NULA(k)

Figure 5 shows DULA [START_REF] Baeza-Yates | A faster algorithm for approximate string matching[END_REF], that has 8 states. DULA(2), the automaton obtained from NULA(2) visible on Figure 4, has 50 states, and is not represented here. For each value of k, the number of states of DULA(k) is computed with R(k).

R k (x, y) = (x ,y )∈Q k ,y<y ,(x,y) (x ,y ),(x ,y ) (x,y) R k (x , y ) R(k) = (x,y)∈Q k R k (x, y) R k (x, y)
is the number of reduced subsets of Q k such that the state (x, y) is the element with the smallest y in the subset. The total number of reduced subsets for Q k , R(k) is obtained by summing over all possible states of Q k . For k ranging from 1 to 10, R equals 8, 50, 322, 2187, 15510, 113633, 853466, 6536381, 50852018, 400763222. We conjecture that R(k) is in O(7 k ).

•11

100 11• 0 1 1 101 1 1 0 0 1 0 •0 0 111 •1• 1•1 •1• ••1 1•• •01 • 1 0 10• 1 • 0 0 • 1 0 1 • •01
Fig. 5. DULA(1), Deterministic Universal Levenshtein Automaton for k = 1. It has 8 states, numbered from 0 to 7, that are all accepting. In the table, we report the subset of Q1 corresponding to each state, as well as the semantics of the state.

Application to the Neighbourhood Counting Problem

We now turn to the problem of computing the cardinality of Lev(P, k). This value depends on the length of the word, the input alphabet and the internal structure of the word. Consider for example the three-letter alphabet {A,B,L}.

Lev( AAA,1) = { AAA, AA, AAB, AAL, ABA, ALA, BAA, LAA, AAAA, BAAA, LAAA, ABAA, ALAA, AABA, AALA, AAAB, AAAL} Lev(LAB,1) = { LAB, LA, AB, LB, AAB, BAB, LBB, LLB, LAA, LAL, ALAB, BLAB, LLAB, LAAB, LBAB, LALB, LABA, LABB, LABL}

In the first case, the neighbourhood has 17 elements, and in the latter case 19 elements. The combinatorics is even more complex for greater values of k.

We show how to solve this problem efficiently with the help of DULA(k). Indeed, it is enough to intersect DULA(k) with the set of all sequences of bit vectors that are a valid encoding for some string V with respect to P . We designate by Encod(P, k) this latter language.

3.1 A DFA for Encod(P, k) Definition 11. Let P be a string of Σ * . Encod(P, k) is the set {u ∈ ({0, 1} 2k+1 ) * ; ∃V ∈ Σ * s.t. u is the k-encoding of V wrt to P }.

From Definition 2, we know that the elements of Encod(P, k) are strings of k + m bit vectors.

Definition 12. Let V ∈ (Σ ∪ {$}) * . Define B(V ) = {χ(s, V ); s ∈ Σ}.
B(V ) is the set of all bit vectors u of length |V | that satisfies the three following properties.

-

If V i = V j , then u i = u j .
-If u i = 1 and u j = 1, then

V i = V j . -If V i = $, then u i = 0.
Example 13. B(ABL)={000, 001, 010, 100}, B(ABB)={000, 011, 100}, B(AA$)= {000, 110}.

Encod(P, k) is recognized by the following DFA.

the input alphabet is ({0, 1} 2k+1 ) * , the set of states is {0, . . . , m + k} ∪ {$ m-k+1 , . . . , $ m+k }, the transition function γ is defined by

γ(i -1, u)=i, 1 ≤ i ≤ m + k, u ∈ B(P [i -k..i + k]) γ(i -1, 0 k+1+m-i 1 i+k-m )=$ i , m -k + 1 ≤ i ≤ m + k γ($ i-1 , 0 k+1+m-i 1 i+k-m )=$ i , m -k + 1 < i ≤ m + k
the start state is 0, there are two accepting states: m + k and $ m+k . Each state i recognizes the encodings of strings of Σ + of length i, and each state $ i recognizes the encodings of strings of Σ + $ + of length i. The transition from i -1 to $ i corresponds to the first occurrence of $ in the string. Figure 6 shows the DFAs obtained for AAA, LAB and k = 1. We also give the DFA for BALLAD and k = 2, for which several encodings were provided in Example 3.

Back to the counting problem

Considering that Lev(P, k) = DULA(k)∩Encod(P, k), we will exploit the product automaton of DULA(k) and Encod(P, k). This product automaton needs not to be constructed explicitely. It will serve us to design the recurrence formula to compute the size of Lev(P, k).

From the product automaton, we can deduce what is the size of the language recognized by Lev(P, k). This is simply the total number of distinct paths leading from the start state to an accepting state. What we still have to do is to bring the problem back to the initial alphabet Σ. For that, we need a function α : {0, 1} 2k+1 × Σ 2k+1 → IN that computes the number of symbols s of Σ such that χ(s, V ) = u, for each bit vector u and each word V over Σ. α(u, V ) = 1, whenever at least one bit of u is 1 α(u, V ) = |Σ| -|V | Σ otherwise (in this case, u = 00 • • • 00) α is used to assign a multiplicity to each transition of the product automaton. In this context, the total number of underlying strings of Σ * is the sum of all multiplicities of all distinct paths. Define S as follows.

S : Reduced(Q k ) × ({0, . . . , m + k} ∪ {$ m-k+1 , . . . , $ m+k }) → IN S(0, 0) = 1 S(q , i + 1) = u,q,q =δ k (q,u) α(u, P [i -k..i + k]) × S(q, i), 0 ≤ i < m + k S(q , $ m-k ) = u,q,q =δ k (q,u) S(q, m -k -1) S(q , $ i+1 ) = u,q,q =δ k (q,u) S(q, i -1) + S(q, $ i ), m -k ≤ i < m + k
m is the length of P . S(q, i) is the number of distinct paths leading from the start site to the state (q, i) in the product automaton. The final result is obtained by summing over the accepting states. Proposition 14. The labels on the transitions are the bit vectors. For each bit vector, we also mention the corresponding symbol(s) of Σ. The number of such symbols is α. The result is the sum of S over the last column.

|Lev(P, k)| = q∈Reduced(Q k ) S(q, m + k) + S(q, $ m+k ) A 1 1 1 1,$ 3 0,$ 4 AAA 1 1 2,$ 4 0,1 6
Figure 7 shows the developments of S for Lev(AAA, 1) and Lev(LAB, 1), which were obtained from DULA(1) in Figure 5 on the one hand, and Encod(AAA,1), Encod(LAB,1) in Figure 6 on the other hand. S can be implemented by dynamic programming with a table of size R(k) × (m + 3k). Each element of the table is computed in constant time. So the algorithm has a time complexity of O(m). In practice, for each possible word structure (e.g. AAA, AAB, ABL), the associated transitions in DULA(k) can be extracted during a preprocessing step. As for the space complexity, this is not necessary to store the full dynamic table. At any instant, the algorithm only requires elements from the current row (i) and the previous row (i -1). Thus the space complexity is in O(1).

Remark 15. In Section 2, we have put a lot of efforts into defining a universal automaton that is able to process any pair of strings. In this Section, we have specialized this automaton for a fixed pattern P and have obtained a DFA for Lev(P, k). Alternatively, we could have directly used the DFA presented in [START_REF] Schulz | Fast string correction with levenshtein automata[END_REF]. The construction of this DFA is also in linear time, for a fixed value of k. However it requires a complex table-based preprocessing, that is dependent on each pattern. In our approach, the computational burden is reported to the design of DULA(k), which is performed only once. This is a new route to recover the result established in [START_REF] Schulz | Fast string correction with levenshtein automata[END_REF].

Conclusion

We have shown how to count the number of strings present in the neighbourhood of some fixed reference word P . The algorithm produces a product automaton, which could also be used to generate the set of all strings in the neighbourhood of P , or to sample it. This generic approach could extend to other target regular languages, instead of the singleton language {P }. The downside of the method, however, is that it requires the computation of the DFA DULA(k), whose size increases exponentially with the number of errors k. It exceeds one million states with k = 8. Another route is to construct the product automaton with NULA(k), instead of DULA(k), and then determinize the resulting automaton with an optimization similar to Proposition 9. This could lead to a lower memory consumption for regular languages that do not need all transitions of DULA(k).

Fig. 2 .

 2 Fig. 2. Bit vector representation for the string SALAD with respect to BALLAD

Example 3 .

 3 2-encodings with respect to the pattern BALLAD. 10010 00010 00011 00111 01111 B A L D 00100 00100 00110 00001 00001 00011 00111 01111 B A L L 00100 00100 00110 00110 00001 00011 00111 01111 B A L L A D S 00100 00100 00110 01100 00100 00100 00000 01111

Definition 5 .

 5 Let k be a positive number. The Nondeterministic Universal Levenshtein Automaton for k, denoted NULA(k), is the NFA represented as follows.

7 00000Fig. 6 .

 76 Fig.6. DFAs for Encod(AAA, 1), Encod(LAB, 1) and Encod(BALLAD, 2)

Fig. 7 .

 7 Fig.7. Computation of S for AAA and LAB, Σ = {A, B, L} and k = 1. For each state (q, i) or (q, $i) (left part of the box), we indicate the value of S (right part of the box). The labels on the transitions are the bit vectors. For each bit vector, we also mention the corresponding symbol(s) of Σ. The number of such symbols is α. The result is the sum of S over the last column.


Anthony Boyer

Hugues Duffau

Emmanuel Mandonnet

Marion A Vincent

Sofiane Ramdani

David Guiraud
email: david.guiraud@inria.fr

François Bonnetblanc

Francois Bonnetblanc
email: francois.bonnetblanc@inria.fr

Attenuation and delay of remote potentials evoked by direct electrical stimulation during brain surgery

Keywords:

Introduction

Direct electrical stimulation (DES) is used to perform intraoperative functional mapping of the brain and guide tumor resection during awake neurosurgery of low-grade gliomas. It generates transient behavioral/cognitive disturbances allowing the identification of both cortical areas and subcortical white matter pathways, which are essential to the function [START_REF] Duffau | Stimulation mapping of white matter tracts to study brain functional connectivity[END_REF][START_REF] Mandonnet | Direct electrical stimulation as an input gate into brain functional networks: principles, advantages and limitations[END_REF]. DES electrophysiological effects remain poorly understood and the recording of Evoked Potentials (EPs) may be used to study the conductive and integrative properties of the neural ensembles being (in)directly affected [START_REF] Vincent | The difference between electrical microstimulation and direct electrical stimulationtowards new opportunities for innovative functional brain mapping?[END_REF]. Three types of cortical EPs should be differentiated: (1) the Direct Cortical Response (DCR), recorded in the immediate vicinity (same gyrus) of the cortical stimulation site and commonly occurring with a 20 ms delay; (2) the Cortico-Cortical Evoked Potential (CCEP), elicited by the physiological propagation through white matter tracts from the stimulated cortical area towards the remote recording site, around 30 ms after the stimulation and (3) the Axono-Cortical Evoked Potential (ACEP), recorded cortically after the stimulation of underlying white matter tracts [START_REF] Vincent | Electrophysiological brain mapping: Basics of recording evoked potentials induced by electrical stimulation and its physiological spreading in the human brain[END_REF].

The DCR typically consists of a negative deflection called the primary negative potential (N1) which peaks between 15 and 25 ms after the stimulation onset and is believed to result from the spatiotemporal summation of post-synaptic potentials at the apical dendrites of pyramidal cells. A stronger stimulus intensity may elicit more complex EPs with spikes of positive polarity and brief duration (approximately 5 ms) appearing in the initial phase of the response, and followed by the N1 component. This early positive deflection, the P0 component, has been ascribed to serial "all or none" discharges in the soma of cortical pyramidal neurons and could be illustrative of the summation of synchronous action potentials [START_REF] Goldring | Direct response of human cerebral cortex[END_REF][START_REF] Li | Cortcial Intracellular Synaptic Potentials and Direct Cortical Stimulation?[END_REF][START_REF] Goldring | Distinctive electrophysiological characteristics of functionally discrete brain areas: a tenable approach to functional localization[END_REF][START_REF] Vincent | Electrophysiological brain mapping: Basics of recording evoked potentials induced by electrical stimulation and its physiological spreading in the human brain[END_REF]). On the other hand, the electrogenesis of CCEP and ACEP remain ambiguous and the "spread of activity in the cerebral cortex" as initially investigated by [START_REF] Bonnetblanc | The spread of activity in the cerebral cortex[END_REF] is mostly unknown [START_REF] Matsumoto | Functional connectivity in the human language system: a cortico-cortical evoked potential study[END_REF][START_REF] Yamao | Intraoperative dorsal language network mapping by using single-pulse electrical stimulation[END_REF].

DES coupled with the recording of EPs could be used to assess the electrophysiological status of cortical sites of interest and more particularly to probe the spatiotemporal connectivity in short-and long-range networks. We hypothesized that the amplitude and latency of EPs may be dependent on the propagation mechanisms (i.e. intracortical vs. subcortical). More specifically, the N1 should be attenuated and, more importantly, delayed, when measured in a different gyrus or lobe because of the additional conduction time and divergence of association fibers.

Material and methods

Two patients underwent an awake brain surgery of a low-grade glioma. The first patient (P1) was a 41-years old man with a glioma located within the left fronto-temporoinsular region. The second patient (P2) was a 50-years old woman with a glioma located within the left paralimbic region. The study (UF 965, n° 2014-A00056-43) was approved by the local ethics committee and patients signed an informed consent. ECoG data were recorded intra-operatively, after the tumor resection and under general anesthesia. Two 4-electrode ECoG strips (2.5 mm platinium contacts, 10 mm spaced, DIXI, France) were positioned on the brain surface of both patients (see B and D). ECoG signals were recorded using differential configurations for both strips, sampled at 10 kHz (PowerLab, ADInstrument) and band-pass filtered from 0.5 Hz to 1 kHz. The front-end amplifier (g.BSamp, G.tec, Austria) was grounded to a patch-electrode located on the right acromion and signals were recorded with a gain of 1,000. Two functional sites identified during the awake mapping were stimulated using biphasic constant-current square wave delivered thanks to a bipolar probe (0.5 mm diameter electrode tips, 5 mm apart) (Nimbus, Innopsys, France). At each stimulation site, 3 stimulation conditions were tested and repeated twice consecutively. The choice of the stimulation parameters was based on previous work which proposed that shortened Pulse Width (PW) and increased Inter-Electrode Distance (IED) may help recruiting different neural elements [START_REF] Boyer | Electrophysiological Activity Evoked by Direct Electrical Stimulation of the Human Brain: Interest of the P0 Component[END_REF]): (1) PW = 1 ms, IED = 0.5 cm; (2) PW = 0.5 ms, IED = 0.5 cm; (3) PW = 0.5 ms, IED = 1 cm. The stimulation intensity was set to 2 mA and frequency was chosen so it allowed long enough time-window to record EPs, especially the N1 component, while keeping the stimulation duration as short as possible (10 and 9 Hz for P1 and P2 respectively). Mean EPs traces were obtained by averaging ECoG signals time-locked to the DES onset. The number of stimuli available depended on the stimulation frequency and its duration, which ranged from 4 to 8 s. The baseline of individual stimulus was set before averaging, using the mean value of the last 10 ms preceding each DES artifact. Also, mean EPs obtained using the first and last 15 repetitions were compared in order to verify for EP stability over time, controlling for possible bias induced by the 9-10Hz stimulation, especially progressive polarization. Finally, 99% confidence intervals (Standard Error of the Mean) were estimated for each DCR as they showed measurable EPs after each stimulus.

Results

The experiment lasted 12 and 9 minutes for P1 and P2 respectively. In total, 12 stimulations were performed on each patients (3 stimulation conditions, repeated twice, on 2 sites) and cortical responses were recorded between each electrodes of both strips (6 channels measuring differential signals) resulting in 72 traces per patient. 35% of traces (31 for P1 and 19 for P2) were not exploitable because of faulty measurement and/or very poor signal-tonoise ratio. 30% of traces (22 for P1 and 22 for P2) did not show noticeable EPs after averaging. Finally, 35% of traces (19 for P1 and 31 for P2), showed EPs after averaging, and DCRs were sometimes distinguishable on raw signal.

Overall, all measured EPs were similar in shape but delayed in time and largely attenuated in amplitude when elicited from a different gyrus or lobe than the one being recorded. More importantly, we were able to measure EPs triggered by both close and remote stimulations for a few recording locations. The recording channels discussed and presented in Fig. 1, are the ones for which we were able to record both DCRs and CCEPs, elicited by intraand extra-gyral DES respectively. This was the case for P1 when stimulating Wernicke's area (S1) which produced local DCRs and elicited CCEPs over the ventral premotor cortex and inversely (see panels A and B). For P2, we recorded DCRs and CCEPs over the precentral gyrus after stimulating respectively the ventral premotor cortex (S2) and the middle part of the superior temporal gyrus (S1), despite the presence of the Sylvian fissure and the operative cavity in-between (see D). Interestingly for this patient, we also measured CCEPs over the superior temporal gyrus when stimulating the precentral gyrus (S2) but abnormally delayed EPs when stimulating residual temporal tissue next to the operative cavity (S1) (see C). Regarding stimulation parameters, we were able to trigger long distance CCEPs when using 1 cm IED with regard to other stimulation conditions which may not induce EPs (see B and D) but they also appeared altered in both latency and amplitude when compared to other existing CCEPs (see A and C).

P1 and P2 brain mappings: Pictures illustrating the stimulation sites (S1, S2) and ECoG positioning with respect to the initial 60 Hz cortical brain mapping (numbered paper tags). Electrodes of both ECoG strips are numbered from 1 to 4 and from 5 to 8. The Sylvian fissure and central sulcus are highlighted by a white dashed lines and annoted "SyF" and "Cs" respectively. For P1, experimental DES was applied on: (1) the Wernicke's area (S1), associated with complete anomia; (2) the ventral premotor cortex (S2), which led to movement and counting interruptions. Strip 1 spans over both temporal and parietal lobe with: electrode 1 over the most posterior part of the superior temporal gyrus; electrode 2 over the Sylvian fissure; electrodes 3 and 4 over the adjacent supramarginal gyrus. Strip 2 spans over the precentral gyrus with: electrodes 5 to 7 over the ventral premotor cortex; electrode 8 is bordering with the most posterior part of the partially resected dorsolateral prefrontal cortex. For P2, experimental DES was applied on: (1) the middle part of the superior temporal gyrus (S1) which led to complete anomia; (2) the precentral gyrus (S2), which induced articulatory disorders. Strip 1 spans over the superior temporal gyrus with: electrodes 1 and 2 over its middle third; electrodes 3 and 4 over its most posterior part. Strip 2 spans over the precentral and dorsolateral prefrontal gyri with: electrodes 5 and 6 over the ventral premotor cortex; electrodes 7 and 8 are respectively bordering and within the adjacent dorsolateral prefrontal cortex. Tumor was about 164 cm 3 for P1 and 150 cm 3 for P2. The number of averaged stimuli is reported within parentheses for each trace. 99% confidence interval estimated for DCRs are represented by grey surfaces to demonstrate that CCEPs do not belong to them. Additional traces corresponding to variations of stimulation parameters were added if available, regardless of the presence of EPs. A: Differential recordings between electrodes 2 and 3 for P1 while stimulating S1 (-170 µV, 21 ms delay) and S2 (amplitudes ranging from -40 µV to -17 µV, delays ranging from 25 ms to 38 ms). EPs following S2 stimulation are CCEPs because of the presence of the central fissure between the stimulation and recording sites. The EP measured after stimulating S1 is ambiguous because electrode 2 lies on the Sylvian fissure, but the short latency and enhanced amplitude with regard to the CCEPs suggest a DCR. Note the dashed line indicating a different amplitude scale for the DCR, which was reduced by a factor 3 for visualization purposes. B: Differential recordings between electrodes 5 and 6 for P1 while stimulating S2 (-75 µV, 20 ms delay) and S1 (-44 µV, 30 ms delay). EP following S1 stimulation is a CCEP because of the presence of the Sylvian fissure between the stimulation and recording sites. EP following S2 stimulation should be viewed as DCR as it was recorded on the same gyrus and it showed shorter latency and enhanced amplitude in comparison with the CCEP. C: Differential recordings between electrodes 3 and 4 for P2 while stimulating S1 (amplitudes ranging from +29 µV to +62 µV, delays ranging from 52 ms to 62 ms) and S2 (+36 µV, 32 ms delay). EP following S2 stimulation is a CCEP because of the presence of the Sylvian fissure between the stimulation and recording sites. EPs following S1 stimulations should be viewed as DCRs as they are recorded on the same gyrus but the latencies and amplitudes appeared unusual. EPs are positive because of differential measure. D: Differential recordings between electrodes 6 and 7 for P2 while stimulating S2 (amplitudes ranging from -260 µV to -310 µV, 20 ms delay) and S1 (-24 µV, 38 ms delay). EP following S1 stimulation is a CCEP because of the presence of the Sylvian fissure and the operative cavity between the stimulation and recording sites. EPs following S2 are likely DCRs as they are recorded on the same gyrus, which is corroborated by their short latencies and maximized amplitudes with regard to the CCEP. Note the dashed lines indicating different amplitude scales for the DCRs, which were reduced by a factor 5 for visualization purposes.

Discussion

Variations in amplitude and delay of EPs are most likely due to different propagation mechanisms, which can be intra-or sub-cortical, and correspond to previously described DCRs and CCEPs. DCRs were not necessarily recorded in the immediate vicinity of the stimulation site (less than 2 mm as stated by [START_REF] Goldring | Direct response of human cerebral cortex[END_REF] and can be observed at greater distances, as long as both the stimulation and recording are performed on the same gyrus. The electrode may thus detect the intracortically propagated response initiated at the stimulation site. By contrast and by definition, CCEP should be observed when DES and recording sites are not in the same gyrus. CCEP recorded in this experiment were delayed and attenuated in comparion with DCRs. The attenuation for CCEPs is probably the result of the limited temporal and spatial summation of post-synaptic potentials (PSPs) consequent to the divergence of white matter pathways arising out of the distant stimulation site. This contrasts with DCR, which involves the synchronous summation of PSPs in the vicinity of DES. Increased delays, on the other hand, may be caused by extra propagation time and be indicative of the conductive properties of the [short] association fibers being recruited. Unusually delayed EPs noticed for P2 when stimulating and recording the superior temporal gyrus may be evocative of pathologically altered fibers (but subcortical intragyral pathways are unknown) or implies intracortical propagation of DCRs which may be delayed and attenuated with distance. These late responses may also be interpreted as further integrative processes. Particular attention should be paid to the delays and amplitudes of EPs when the recording is performed over adjacent gyrus, because of the possible overlap of short latency CCEPs, propagated along intergyral association U-fibers, and DCR, elicited because of local electric conduction of the tissue. In such situation, it could be difficult to disentangle DCRs from short CCEPs. However, there are no such ambiguities with the EPs we presented as the stimulation and recording sites are on different lobes when measuring CCEPs. As a whole, differentiating CCEP from DCR based on anatomical landmarks is problematic. Due to a limited spatial resolution of most recording systems and methods, it is somehow difficult to identify the exact generator mechanisms of the evoked responses that may combine together. In consequence, the nomenclature for responses adjacent to the DES site is not yet totally clear. Indeed, they are also sometimes called "adjacent large CCEP responses" [START_REF] Shimada | Impact of volume-conducted potential in interpretation of cortico-cortical evoked potential: Detailed analysis of high-resolution electrocorticography using two mathematical approaches[END_REF][START_REF] Kobayashi | High frequency activity overriding cortico-cortical evoked potentials reflects altered excitability in the human epileptic focus[END_REF] while some prefer the classic term DCR.

All EPs we measured were very similar in shape, suggesting stereotyped electrophysiological response triggered by massive and synchronous activation rather than genuine output integration. The summation of PSPs is the key element shaping ECoG signal and appears loosely dependent of that occurring in input. According to previous works focusing on the peripheral system, the IED and PW were tuned in order for the electrical field to reach deeper and bigger neural elements in the cortical column and facilitate the generation of action potentials with less charge injections, namely to improve the generation of CCEP [START_REF] Vincent | Electrophysiological brain mapping: Basics of recording evoked potentials induced by electrical stimulation and its physiological spreading in the human brain[END_REF]. We observed a few differences depending on DES parameters, especially when using 1 cm IED, which produced more consistent and earlier, but altered in amplitude, CCEPs. This is in adequacy with some preliminary results showing that increased IED may elicit the P0 component of the DCR, suggesting a greater spatio-temporal summation of action potentials propagating on sub-cortical association fibers [START_REF] Boyer | Electrophysiological Activity Evoked by Direct Electrical Stimulation of the Human Brain: Interest of the P0 Component[END_REF].

The use of a 9-10 Hz stimulation may prevent the recording of later components which are observed in more recent CCEP studies performed on epileptic patients and using 1 Hz stimulation [START_REF] Matsumoto | Functional connectivity in the human language system: a cortico-cortical evoked potential study[END_REF][START_REF] Yamao | Intraoperative dorsal language network mapping by using single-pulse electrical stimulation[END_REF]). However and importantly, there are major differences in DES parameters between these studies and older ones (which focused on DCR), and ours [START_REF] Vincent | The difference between electrical microstimulation and direct electrical stimulationtowards new opportunities for innovative functional brain mapping?[END_REF][START_REF] Vincent | Electrophysiological brain mapping: Basics of recording evoked potentials induced by electrical stimulation and its physiological spreading in the human brain[END_REF]. This is especially true when considering the current intensity. In our study we wanted to investigate whether classical N1 components could be observed with lower intensity of DES in order to use these electrophysiological responses in real-time, as putative indicators of connectivity and perform electrophysiological mapping during awake brain surgery. Interestingly, [START_REF] Goldring | Direct response of human cerebral cortex[END_REF] demonstrated that the shape of DCR was not altered when DES frequency was inferior to 20Hz, but a slownegativity appeared (whose frequency was inferior to 1 Hz). As for our data, we verified the stability of EPs over the stimulation period comparing subsets of averages.

From a methodological perspective, differential configuration allowed the recording of signals with better signal-to-noise ratio and more localized detection, at the expense of an unconventional shape of EPs. It remains particularly difficult to control for the electronic noise in such a neurosurgical context and a large part of the data has to be discarded. Nonetheless, it seems possible to discriminate DCRs from CCEPs and explore anatomofunctional connectivity using such configuration. We were, for instance, able to observe the bidirectional nature of the arcuate fasciculus as described in [START_REF] Matsumoto | Functional connectivity in the human language system: a cortico-cortical evoked potential study[END_REF].

It is important to note that these measurements were performed on patients under general anesthesia, at the end of the surgery. Anesthesia may play some role on the excitability level and on the shape of the response. More specifically, it has been shown that drugs may be correlated with an increase of the N1 component and the disappearance of the (rare) N2 component [START_REF] Goldring | Direct response of human cerebral cortex[END_REF]. However, the latencies and the amplitude of the N1 responses were not modified consistently by drugs injection.

Conflict of interest statement

On behalf of all authors, the corresponding author states that there is no conflict of interest.

Sponsorships and funding

Research supported by the LabEx NUMEV project (n° ANR-10-LABX-20) funded by the French government's "Investissements d'Avenir" program managed by the French National Research Agency (ANR) and complementary grants from the Institut Universitaire de France and INSERM laboratory (U1093).


PHD Shasha Shi
email: shasha.shi@telecom-bretagne.euno

Patrick Gioia
email: patrick.gioia@orange-ftgroup.com

Gérard Madec
email: gerard.madec@telecom-bretagne.euno

High efficient computational integral imaging reconstruction based on

Keywords: integral imaging, computational reconstruction, parallel-group projection, sub image

come

Introduction

3D technology has received considerable attention since 1838 when Charles Wheatstone first described his theories on stereopsis. In spite of its significant progress for past decades of years, 3D technology has not yet accomplished its promise of offering natural viewing effect for spectators. Currently provided stereoscopic or auto-stereoscopic systems produce relief perception to viewers by simulating binocular disparity between two eyes. In these systems, slightly different images are sent to left eye and right eye separately to construct a 3D perception for viewers. Due to its relative simplicity, this kind of 3D system is widely utilized for today's 3D products but it contains some inevitable defects such as the discontinuous viewpoint for different viewing directions and irrespective effect between accommodation and convergence operations [START_REF] Maritinez-Corral | Progresses in 3D Integral Imaging with Optical Process[END_REF].

According to the inborn limitations of conventional stereoscopic and auto-stereoscopic systems, next generation 3D systems providing parallax in all directions, make accommodation and convergence work in unison, are eagerly required. Integral imaging is proposed as a practical and prospective technique for this purpose.

A typical integral imaging system could be dived into two stages, named as acquisition and display. Figure 1-1 illustrates the optical working principle for the integral imaging system. In the pickup process, the direction and intensity information of light rays come from 3D object are captured through an array of convex lenses or pinholes and then recorded by a 2D image sensor, as shown in the left picture of Figure 1-1. The recorded image through each lens or pinhole is so called elemental image (EI), which represent the 3D objet from a specific perspective. An integral image (II) is composed of many EIs, with their number corresponding to the number of lenses or pinholes. While in the display side, the decoded array of elemental images is exposed behind the same optical devices (lenses or pinholes array) to reconstruct the light rays of the 3D objet in the reverse direction, shown in the right picture of Figure 1-1. For optical integral imaging reconstruction (OIIR) operation, the diffraction and restrict of optical devices introduced the quality degradation and the viewing angle limitation of reconstructed image. To overcome these problems, computational integral imaging reconstruction (CIIR) method has been developed recently [2, 6 and 7], the CIIR is achieved by simulating the projection effect of optical rays of EIs in OIIR using a computer built virtual pinhole. In this method, the reconstructed output image (ROI) could be generated at any arbitrary distance from virtual pinhole array without device degradation and optical diffraction effect. Although computational integral imaging system contains several useful superiorities, there still exits some limitations which obstruct its practical utilization, such as low ROI quality and high calculation complexity. In this paper, a new CIIR approach based on parallel-group projection (PGP) method was proposed to solve these problems.

This paper is organized as follows. Section 2 presents some related works on computational integral imaging reconstruction. The proposed CIIR method based on Parallel-Group Projection is explained in section 3. Simulation results and detailed analysis are reported in section 4. Finally, conclusions and discussion are drawn in section 5.

Related work

Figure 2-1 illustrates the process of conventional CIIR method. In this method, each elemental image is inversely projected through its corresponding pinhole, magnified by the magnification factor M, which is the quotient between the distance from the virtual pinhole array to the reconstructed image plane (denoted by 'z' in figure 2-1) and the distance from elemental image plane to pinhole array (denoted by 'g' in figure 2-1). The magnified pixels are overlapped by each other in reconstructed image plane to reproduce the ROI in the determined distance [START_REF] Seung-Hyunhong | Three-dimensional volumetric object reconstruction using computational integral imaging[END_REF].

Series of reconstructed image with different 'z' value are generated continuously to compose a volumetric reconstruction of the 3D scene. This CIIR approach is straightforward but it contains a crucial problem about image quality degradation. This quality degradation comes from the magnification process, where each pixel is enlarge and superposed by another pixel, and the superimposition between adjacent pixels causes the undesirable interference result in output image.

Figure 2-1 Reconstruction principle of traditional CIIR method

To solve this quality degradation problem, another CIIR method using pixel to pixel mapping technique was proposed in reference [START_REF] Shin | Computational integral imaging reconstruction method of 3D images using pixel-to-pixel mapping and image interpolation[END_REF], the principle of this approach is described in figure 2-2. It removes the magnification process of previous CIIR techniques, generates a dot-pattern ROI by projecting each pixel of the elemental images into output plane, and then fulfills the empty pixels in this dot-pattern image by the linear interpolation method. Because of the exact generation of pixel position, this method can remove the interference effect betweens adjacent pixels and improve the quality of reconstructed image. However, due to the repetitive projecting of each pixel in elemental images, the calculation complexity is still very important. Moreover, as a result of the separated pixel projection process, this method induced an undesirable dotdiscrete effect for ROI. To overcome all the inconvenient of previous techniques, we proposed a new CIIR approach using Parallel-Group Projection (PGP) method in this paper. Different from exited techniques, this new approach propagates sub image (SI) instead of original elemental images through the computer synthesized pinhole array to reconstruct the ROI, where SI is a special form of 2D image which records the information of parallel light rays from the given 3D scene [START_REF] Park | Three-dimensional Optical Correlation Using a Sub-image Array[END_REF], so it is nature that the reverse mapping rays of each pixel from the SI into the reconstructed image are also mutually parallel. Therefore, we considered reconstructing the ROI by parallel-group projection instead of pixel-to-pixel projection to exploit this mapping regularity. This new CIIR technique is able to eliminate the interference and dot-discrete effect in conventional generated ROI and increase the efficiency of previous CIIR methods significantly.

Process of Proposed CIIR Method

A typical scheme of integral imaging system using the proposed CIIR method is illustrated in Figure 3-1. This system can be achieved by three steps. Firstly, sub images recording the light information of parallel rays are generated by OpenSceneGraph (OSG) 3D graphics toolkit using special configured cameras [START_REF] Lee | Design and implementation of a Fast Integral Image Rendering Method[END_REF]. Thereafter, the obtained SIs are encoded and decoded by proper compression technique. Finally, in the display side the output 3D image is reconstructed by proposed CIIR technique. In the following part, we will concern on the concrete procedure of proposed CIIR. As figure 3-2 illustrated, the projected image of the original EIs through a computer synthesized pinhole array could be represented as the accumulation of the mapped result of each SI through the same virtual pinhole array. Due to the known parallel character of the mapping light rays [START_REF] Fellner | Toward the Light Field Display: Autostereoscopic Rendering via a Cluster of Projectors[END_REF], all the mapped pixels from one SI are arranged as a grid with a given interval in the reconstructed plane. Therefore, from the location of one projected pixel we are able to estimate the projected coordinate of every pixel by simply displacing the first mapped pixel with an specific offset, and the eventually ROI is formed by accumulating the projected result of every SI on the reconstructed plane. (1)

ܯ ൌ ሺܼ כ tan ߠ כ 2ሻ/ሺሺ݊ െ 1ሻ כ ሻ (2)
Where 'X ', 'Y 'express the mapped pixel coordinates in ROI. 'ݔ' and 'ݕ' represent the index of the projecting SI, which are ranged as (0, േ1, േ2 … േ ሺn െ 1ሻ/2) . '݊' is the number of SI in horizontal or vertical axis. 'M' is the magnification factor. '' is the size of pixel in SI. 'ߠ' is the viewing angle of the marginal SI Once we get the correct projected position of the first pixel in each SI from formula ( 1) and ( 2), the second step is to deduce the mapped coordinates of the other pixels in this SI through the regular layout of projected pixels, which could be represented simply by the fowling equation:

ሺܺ, ܻሻ ൌ ሺ ܺ , ݍ ܻ ሻ (3)
Where 'p' and 'q' are the position of each pixel in the SI Finally, the projected outcomes of each sub image are superimposed to each other for constructing the displayed ROI at the determined view distance, and a full 3D volume scene is generated by iteratively repeating this process for variant visual distance. Besides its high computational efficiency derived from parallel-group-projection, this CIIR approach contains several other merits compare to previous CIIR methods. First of all, because of its direct application of rendered SI, the transformation and reverse transformation process between II and SIs are no longer necessary in the novel system. Another supplementary advantage lies in the variability of the reconstructed image size, which is in direct proportion to the size of SI, and the size of SI could be adjusted easily by exist interpolation techniques.

Experiments and Results

The experiment was carried out to demonstrate the effectiveness of proposed approach. In the beginning of the experiment, nine sub images with distinguish viewing angles are rendered by OSG toolkit with orthotropic projective cameras, which are illustrated in figure 4-1. The resolution of each SI is 512x512, viewing distance of the virtual camera is 15, number of SI in process with a certain magnification factor, the proposed PGP CIIR method shows a significant efficiency progress compare to traditional and pixel-to-pixel mapping method, and its calculation time is not affect by the increasing of magnification factor M.

Conclusions

In this paper we propose a novel CIIR technique based on parallel-group projection method to improve the efficiency and performance of traditional CIIR methods. In this approach, the reconstructed output image is generated by using the parallel-group projection technique to each sub image. Where sub image is an alternative format of 2D image which recodes the information of parallel light rays from the 3D scene. Exactly and completely pixel position calculation process could avoid the interference problem and the dot-pattern effect of previous CIIR techniques, and the direct projection of SI and the simple mapping regularity of each SI reduce the consuming time of earlier approaches. The experimental results show that our proposed method outperforms the previous CIIR method both in terms of visual quality and in terms of calculation efficiency. In the future, we intend to extend this reconstruction idea from digital visualization system into real optical display system so as to make the present 3D view similar to real life 3D version.

Figure 1 - 1

 11 Figure 1-1 Optical setup for generating (left) and displaying (right) 3D scene in integral imaging system

Figure 2 - 2

 22 Figure 2-2 Reconstruction principle of pixel-to-pixel mapping CIIR method

Figure 3 - 1 Figure 3 - 2

 3132 Figure 3-1 Block diagram of proposed integral imaging system

Figure 3 - 3

 33 Figure3-3 illustrates the detailed calculation process of the proposed CIIR method, which is carried out by three steps. In the beginning, we project the first pixel of a SI into the reconstructed image plane through the computer synthesized pinhole array based on the geometric optical diagram shown in figure3-4. According to this diagram, the projected coordinate of one pixel from SI could be calculated by the following equation: ሺܺ , ܻ ሻ ൌ ሺെݔ כ ,ܯ െݕ כ ܯሻ

Figure 3 - 3 Figure 3 - 4

 3334 Figure 3-3 Block diagram of proposed CIIR method

  axis is 3, and the viewing angle of marginal SI is 15°. According to the known parameters and formula (2), we can get the magnification factor M equal to 25.

Figure 4 Figure 4 - 2 Figure 4 - 3

 44243 Figure 4-1 Generated SIs for a 3D scene

Table 4 -

 4 1 Comparison of the computation cost of three CIIR methods

	Magnification	Traditional CIIR	CIIR using pixel-to-	CIIR using parallel
	parameter		pixel mapping	group projection
	15	67078	703	328
	35	363218	703	328
	55	897578	797	328
	75	1688047	859	328
	95	2724593	922	328
	115	4004547	968	328
	Three different reconstructed images are generated separately by traditional, pixel-to-pixel
	based and parallel-group projection based CIIR methods. The performance of CIIR technique is
	evaluated in terms of image quality and calculation efficiency. Figure 4-2 shows the reproduced
	ROI of three techniques with identical parameters. Compare to traditional methods, proposed
	CIIR approach provide a high quality reconstructed 3D image without interference effect in
	traditional method and point discrete effect in pixel-to-pixel mapping method.

Table 4 -

 4

1 and figure 4-3 illustrate the computational cost of the three methods in function of the magnification factor M. According to the exposed data, we can detect that for a reconstructing


Rémi Bertrand
email: remi.bertrand@lcts.u-bordeaux.fr

Gerald Camus
email: camus@lcts.u-bordeaux.fr

Olivier Caty
email: caty@lcts.u-bordeaux.fr

Vincent Mazard

Sébastien Denneulin
email: sebastien.denneulin@safrangroup.com

Patrick Weisbecker

Jérôme Pailhès

Francis Rebillat
email: rebillat@lcts.u-bordeaux.fr

Vincent Mazars
email: mazars@lcts.u-bordeaux.fr

Jérôme Pailhes
email: jerome.pailhes@ensam.eu

In-situ

Keywords: Tomography, in-situ tensile tests, oxidation, Ceramic matrix composites, Self-healing behaviour

Introduction

The field of civil aeronautics seeks for replacing some of the structural components made of metallic alloys present in the engines [START_REF] Naslain | Boron-bearing species in ceramic matrix composites for long-term aerospace applications[END_REF][START_REF] Habarou | Applications of CMCs to Turbojet Engines: Overview of the SEP Experience[END_REF]. During their use, these parts undergo high temperatures and relatively high stresses, under oxidizing and corrosive environments. SiC/SiBC composite with their fair resistance to high temperatures and excellent mechanical properties are candidates for such applications. However, when the applied stress is sufficiently high, a crack network is created within the material. Each of these microcracks is paths of oxygen to the carbon interphase and the SiC fibre. Oxidation of the interphase (which is essential to avoid a brittle behaviour) and of the fibre drastically reduces the lifetime of the material since 400°C [START_REF] Filipuzzi | Oxidation mechanisms and kinetics of SiC-matrix composites and their constituents[END_REF]. Different types of matrices have been developed to limit oxidation of the interphase [4-5-6-7]. A matrix that combines carbides (deposited by CVI) has thus been especially developed with specific sequences of Si, C, B which generate protective oxides for a wide range of temperature [START_REF] Buchanan | Permittivity and AC conductivity in Yttria-stabilised zirconia[END_REF][START_REF] Baroumes | An Improved Long Life Duration CMC for Jet Aircraft Engine Applications[END_REF]. In particular, boron oxide exhibits a low melting point (almost 460°C), and an appropriate viscosity which can fill up matrix cracks at low temperature (almost 600°C) [START_REF] Buchanan | Glass sealants for carbon-carbon composites[END_REF]. This matrix can be used to protect the material from room temperature at high temperature because as soon as 450°C the layers with boron oxidize create self-healing glass.

Previous studies have showed that this selfhealing glass can modify the behaviour of the material according to the thermomechanical cycle endure by the composite [START_REF] Bertrand | Influence of thermomechanical cycling on the behaviour of a selfhealing Ceramic Matrix Composite in an oxidizing atmosphere" under review[END_REF]. Indeed, the solid oxide presents in the crack involves a stiffening phenomenon and can, in some cases, prevent the reclosure of the crack during the unloading. Micro computed tomography (μCT) and tensile tests under SEM give the possibility of 3D image and analyse in-situ the material at different resolutions while loading it mechanically and thermally. These analyses give access to a very precise insight of behaviour of the oxide in the cracks [START_REF] Bale | Real-time quantitative imaging of failure events in materials under load at temperatures above 1600[END_REF][START_REF] Caty | Application of X-Ray computed micro-tomography to the study of damage and oxidation kinetics of thermostructural composites[END_REF]. This technique is used here for the in-situ study of damage at different stress states and at different temperatures up to 800 °C on small samples.

Material and methods

2.1) Material

The SEM tensile tests are performed on a microtest tensile stage (Deben Microtest). This tensile stage is has a 1kN load cell and tensile speeds are in the range of 0.02 mm.min -1 to 0.4 mm.min -1 . The material tested is a mini-composite composed of about 500 unidirectional SiC fibres (Hi-Nicalon fibres). In order to improve the protection of the SiC-Hi-Nicalon fibres, 50 µm of matrix of [Si-B-C] is deposed by chemical vapour deposition.

A special tension device (Fig. 1) was developed for in-situ tension tests in μCT scanners. This device is adapted for the GE-VTomx μCT and most of the synchrotron beam lines. The design was inspired from the one developed at INSA de Lyon [START_REF] Salvo | X-ray micro-tomography an attractive characterisation technique in materials science[END_REF] adapted for small μCT scanners using a hydraulic actuator and a joule effect heater. The device is designed for small specimen with a section up to 5x5mm 2 and an utile length of 30 mm max. In this paper the utile length of the sample is 2.5x2.5x10 mm 3 to fully image this section at a resolution of 7 μm. The maximum load is 5000 N and the maximum temperature available with this specimen and the power supply used (16V/16A) is 1250°C. The device is cooled by air, the temperature is measured by a pyrometer and the temperature is stable (+/-1%) in the utile length of the sample.

2.2) SEM and μCT tensile tests

Before SEM tensile test, each sample is damaged and oxidized in an ageing furnace (Fig. 2). This method allows to control the cooling and the unloading. After this, the sample is placed on the tensile stage in the SEM. The tensile test is performed at room temperature and environmental SEM atmosphere. During the test, the evolution of a crack is observed by SEM to follow the solid oxide in the crack when the load is applied. Extensions were designed to be plugged in the ageing furnace and to hold SEM extensions for tensile tests (Fig. 2). These extensions are unmounted by removing the pin.

Tests were designed to analyse healing and in particular crack blockages due to the presence or not of solid oxide in the crack (depending on the type of cooling after or before the unloading). The test sequence is described on the Fig. 3. Initially, the sample is scanned without load at room temperature (scan I). Then the load is applied to create microcracks in the material (scan II). While the load is maintained, the temperature is increased to oxidize and heal the material (800°C). After the cooling under load, the material is imaged maintaining the load (scan III). Before scan IV, load is removed to analyse crack blockages. Cracks are reopened (scan V) and the temperature is increased to 800°C to liquefy oxide. Then, the load is removed before temperature to analyse the reclosing of the microcracks when the oxide is liquid (scan VI). Finally, scan VII is performed with load to observe effect of this expulsion on failure morphology (scan VII). Stress was 150 MPa to damage the material without global failure. The temperature is raised to 800°C to liquefy the oxide without reaching the creep temperature of the SiC phases.

3.1) SEM tensile tests

After the ageing of the sample, it is possible to decrease the temperature before the load. In this case, the solid oxide present in the crack blocks the closure of the crack. When the load is applied under SEM, a failure appears between a side of the crack and the solid healing oxide (highlighted by arrows on figure 3-A).

In the other case, if the cooling is done after the unloading, the liquid healing oxide is ejected out of the microcrack. This is highlighted by an arrow on figure 3.B. The liquid oxide at 800°C has a perfect wettability and it can coat the composite as soon as it is ejected out of the crack. When the load is again applied, some microcracks appear in the solid oxide layer present in the crack (see arrows on the Fig. 3.C).

3.2) μCT tensile test

The Fig. 4-A presents the evolution of a crack in the sample for 7 loading steps. Initially (step 1, 0N), no crack is present in the material. Scan 2 is acquired after loading the material without temperature. A crack appears highlighted by the arrow. This crack evolves during the oxidation phase (step 3). After the first thermomechanical cycle (i.e. cooling under load), the crack is always opened whereas there is no load during the scan (step 4). Illustrating the blockage of the crack by the healing oxide at room temperature. When the load is applied again no difference are visible (step 5). On the contrary, when the load is removed before cooling, the crack is almost totally closed oxide not visible at this resolution on µCT (step 6). The crack is reopened when the load is applied again (step 7). Illustrating the blockage of the oxide as seen on figure 3. Fig. 4-B presents a 3D reconstruction of the material and the propagation of the crack already presented in the figure 6. The green segmentation is the crack propagation during the step 5 (i.e. when the load is maintained). On the contrary, the blue segmentation is the crack during the step 6 (i.e. when the load is removed before cooling). The green segmentation is more developed than the blue. This is consistent with the better reclosure of the crack when the material is unloading before cooling (step 6). The volume of the crack is given in the table Fig 4-C. The volume of the crack to step 4 is almost the same as step 3 and 5. This is consistent with the fact that the material is unloading. However, this volume is still larger than that obtained in step 6.

Discussion

The SEM analyzes were conducted for their accuracy (magnification factor x6000, Fig. 5-B) although on a localized area of a material model (mini-composite). The study by μCT gave comprehensive information on real composites at lower resolutions (7μm). In this paper, to obtain more information on the thermo-mechanical behaviour of the material, SEM and µCT were combined and run with in-situ tensile tests.

Thanks to tensile tests under SEM, it is possible to observe that, according to the type of cycle (cooling after or before the unloading), the material has not the same behaviour when a load is applied again:  The first behaviour is when the cooling occurs before the unloading. As soon as a new load is applied the new crack occurs between solid oxide layer and the edge of the crack (Fig. 3-A). No new microcracks appear in the solid oxide. The interfacial bond between the oxide and the matrix is a weakness element. There is no acoustic emission and stiffening is lower which is explained by an oxide interface / low crack [START_REF] Bertrand | Influence of thermomechanical cycling on the behaviour of a selfhealing Ceramic Matrix Composite in an oxidizing atmosphere" under review[END_REF].

 In the other hand, when the unloading occurs before the cooling, as seen Fig. 3-B, a part of the liquid healing oxide is ejected out of the crack. Then, after cooling, if the load is applied, multicracking occurs in the oxide layer (Fig. 3-C). In this case, the interfacial bond between matrix and solid oxide seems more important that in the case where the cooling occurs under load. To reopen the crack, it is necessary to break the thin solid oxide layer.

The µCT tensile test, analyses in 3D the influence of healing solid oxide on the opening/closing of the crack when the material is unloading throughout composite sample. The solid healing oxide in the crack blocks the reclosure of the crack when the material is cooling under load (step 4 Fig. 4-A). After the inverse cycle, i.e. unloading at high temperature, liquid healing oxide cannot prevent the reclosure of the crack. There is a thin layer of oxide which cracks when the load is applied again. This is supported by the previous tests that show the cracking of the oxide layer (Fig. 3-C). This behavior is accompanied by a stiffening of the material with an emission acoustic which reflects a multiple cracking oxide (Fig. 3-C). This phenomenon is studied in [START_REF] Bertrand | Influence of thermomechanical cycling on the behaviour of a selfhealing Ceramic Matrix Composite in an oxidizing atmosphere" under review[END_REF]. These different in-situ methods of characterization allow us to understand better the behavior of material. The observations and results are qualitative analyse although it is possible to quantify the crack opening during the tensile under SEM (≈ 1.8 ± 0.1µm on the Fig 3-B). The volume of the crack quantified by post treatment segmentation highlights the phenomenon of blockage. Indeed, without load, the crack during step 4 as almost the same volume that step 3 or 5 with load. Moreover, the reclosure of the crack is better when the oxide is liquid because the volume of the crack is smaller than in the other case (step 4, Fig. 4-C).

Conclusion

SEM tensile test and 3D μCT images are used to analyse and visualize crack initiation and propagation in self healing 3D CMC. These analyses were made in-situ using a high temperature tension device specially developed for this material and for SEM and μCT. Two types of tests were carried out. The first allows to observe the influence of the type of cooling (before or after the unloading) on the mechanical behaviour when the solid oxide is present in the microcrack. Depending to the type of cooling, the solid oxide does not crack the same way: (i) if the stress is released after solidification: blocking of the crack, and there is failure along the interface healing oxide / material during the next reloading, (ii) if the material is unloading before solidification: crack closes, and the oxide is ejected. During the next reloading, the rest of oxides will crack and the mechanical behaviour will be degraded. The second test allows to bring to light the behaviour of healing phases submitted to load and temperature. Images presented were treated to better analyse cracks and phases. The cracks are rather straight and perpendicular of the load but located on a side of the sample. At high temperature the healing oxide were also observed and can prevent the reclosure of the crack when the material is cooling under load or be ejected in the other case.

Figure 1 :

 1 Figure 1 : (A) Schematic representation of the ageing furnace (B) Extension for the ageing furnace and the SEM tensile test, (C) Picture of the tension device in the GE/VTomx μCT with a specimen heated at 1000°C

Figure 2 :

 2 Figure 2 : Schematic description of in-situ tension tests at high temperature carried out to study the damage and healing of the material

Figure 3 :

 3 Figure 3 : (A) Illustration of the failure of the healing oxide phase, after cooling under load and tensile test under SEM (B) Liquid oxide partially ejected of the microcrack during the unloading at high temperature (C) microcracks in the solid oxide layer after loading, SEM micrographs

Figure 4 :

 4 Figure 4 : (A) Seven states of the thermomechanical sequence describe in figure 3. A zone of the working portion of the sample. The loading direction is indicated by hold arrow and σ. All acquisition was done at 20°C (B) Two crack segmentations (green after step 5 and blue after step 6). 3D reconstruction of the material after X-Ray computed micro-tomography (C) Volume of the crack as function of the step.

Acknowledgements

This work has been supported by the Centre National de la Recherche Scientifique (CNRS) and Safran Ceramics (SAFRAN GROUP). The authors are grateful to Safran Ceramics for supplying the samples and to Felicien Auguin, Matias Gosset, Bruno Vicin and Jonathan Jouveau for their work during an internship on the development of the device and the in-situ tensile test and the I2M laboratory for the development of the device in particular and Fabrice Prince and Jean Marie Diaz.


Nicolas Bus

Ana Roxin

Guillaume Picinbono

Muhammad Fahad

Towards French Smart Building Code: Compliance Checking Based on Semantic Rules

Keywords: OWL, SWRL, ifcOWL, automated compliance checking, semantic rules, construction regulations, BIM, IFC

HAL is

Introduction

Manually checking models for building code-compliance is time-consuming for architects and construction engineers. According to McGraw-Hill [START_REF]Industry Foundation Classes (IFC) for data sharing in the construction and facility management industries[END_REF], they can spend up to 100 hours on this aspect for complex projects. We, therefore, need algorithms that process construction project information and report non-compliant elements.

Building and construction regulation in France are published as human readable texts.

The translation from texts into algorithms is a complex task that requires construction experts to collaborate with computer scientists. Indeed, texts are often ambiguous, incomplete and sometimes contradictory and contain implicit knowledge. Thus it is crucial to correctly parse and prepare such texts, in order to make sure that all implicit knowledge has been made explicit and that all contradictory statements have been correctly interpreted. CSTB (Centre Scientifique et Technique du Bâtiment, French Research Center for Architecture and Construction) started working on this topic way back in 2002 [START_REF] Maïssa | Regulation checking in a Virtual Building[END_REF]. This resulted in a first mapping between the regulations vocabulary and the IFC vocabulary [START_REF]IFC 2x3, Specification[END_REF].

Today ISO 16739 is the standard for data sharing in the construction and facility management industries. ISO 16739:2013 or IFC (Industry Foundation Classes) specifies a conceptual data schema and an exchange file format for Building Information Model (BIM) data. The conceptual schema is defined in EXPRESS data specification language [START_REF]Industry Foundation Classes (IFC) for data sharing in the construction and facility management industries[END_REF]. While IFC is the standard for exchanging BIM models, other standards exist for describing building systems or elements. Thus checking a building model against compliance constraints has to take into account the other ways in which building data can be described. Given the level of schema heterogeneity of these standards, semantic Web technologies appear as a solution for delivering semantic interoperability among the diversity of schemas. This has resulted in the publication, by bSI (buildingSmart International) of an OWL ontology version of the IFC standard [START_REF]IFC 2x3, Specification[END_REF].

In this paper, we thus present an approach based on semantic Web technologies. We argue that it can bring several benefits when applied to building automated regulation compliance checking.

Related work

Semantic Web Technologies in the context of BIM

In the context of BIM resources modeling has been identified as an interesting approach for achieving information interoperability [START_REF] Mendes De Farias | Semantic Web Technologies for Implementing Cost Effective and Interoperable Building Information Modeling[END_REF]. Existing IFC-related ontologies were conceived as direct syntax mappings between EXPRESS and OWL languages [20][21]. One of the latest and the most solid implementations of an IFC ontology is ifcOWL approach proposed in (W3C Linked Building Data Community Group, 2014). This version of ifcOWL is also a Candidate Standard [START_REF]ifcOWL[END_REF] for build-ingSMART (meaning it is considered as an activity in the process of acquiring international consensus before being submitted to the Standards Committee for a final vote).

Regulation's checking

In the context of the buildingSMART community, the common approach used for verifying IFC files is called Model View Definition (MVD). An MVD "defines a subset of the IFC schema that is needed to satisfy one or many Exchange Require-ments of the AEC industry" [START_REF]Architecture and Construction Workshop[END_REF]. MVDs are serialized in a format called mvdXML [START_REF]MVDXML specification and support implementation[END_REF].

The traditional approach of verifying IFC models using MVDXML is limited as stated in [START_REF] Zhang | Model view checking: automated validation for IFC building models[END_REF] and [START_REF] Mendes | A rule-based methodology to extract building model views[END_REF]. Major identified limitations are its limited scope when applying conditions and constraints on several branches of an IFC model, poor geometric analysis of an IFC model, lack of mathematical calculations, and support of only static verification of a model. Given these limitations, several works as [START_REF] Zhang | Model view checking: automated validation for IFC building models[END_REF] and [START_REF] Mendes | A rule-based methodology to extract building model views[END_REF] investigated approaches based on logical rules as viable alternatives. Pauwels and al.

showed the benefits of using web semantic through a simplified ontology derived from ifcOWL to check construction rules [START_REF] Pauwels | SimpleBIM : from full ifcOWL graphs to simplified building graphs, Ework and ebusiness in architecture, engineering and construction[END_REF]. Fahad et al. studied that building code compliance based on the semantic web rule language has more advantages and cover more functionalities of verification as compared to MVDXML approach [START_REF] Fahad | Towards Validation of IFC Models with IfcDoc and SWRL A Comparative Study[END_REF].

Method

Addressing schema heterogeneity

Given the diversity of schemas that can be defined for containing building data, we started defining a modeling charter. In our vision, a modeling charter is a document describing how construction project stakeholders agree to model and share data among them. This document describes data exchange processes, modeling best practices, export parameters, nomenclatures and classifications to be used. With all actors of the project sharing the same modeling rules drastically decreases the number of possible schemas. The modeling charter is also used to extend IFC properties with adhoc property sets dedicated to specific topics -such as certification level, equipment performances, product impact on environment, etc.

As a future goal, the CSTB aims at providing, in the short term, a national modeling charter core to be used as a common starting point when engineering the modeling plan of a project. We are working in close cooperation with architects, building owner, engineers, audit experts, building quality control agencies and government to reach this goal.

Aligning the Building model with the Regulations ontology

The semantic approach allows describing a building as well as a requirement by using the same atomic fact formalism called triple. In the context of the Semantic Web, a triple represents a statement in the form of "Subject -Predicate -Object", where each part is identified by an URI (Uniform Resource Identifier). Only the object of a triple is an exception to this rule, as it can hold a literal value, and in this case the predicate is called a datatype property. In the case all elements of a triple are identified by URIs, the predicate is called an object property.

The ifcOWL vocabulary [START_REF]ifcOWL[END_REF] includes thousands of classes and properties, very close to building physical components. This vocabulary is not designed to deal with performances, requirements and building functionalities. Thus, if we seek to map this vocabulary to a Regulations ontology (described in section 3.3), we can oversee that the alignment between these two vocabularies cannot be based on a bijective transfer function. Indeed, it is impossible to define a one-to-one mapping between the considered terms. Moreover, some regulatory terms can only be defined with a logical statement aggregating various IFC terms. For instance, in our regulation ontology, a simple regulatory term like "highest storey" is defined by ifcOWL concepts: "IfcBuildingStorey" (ifcOWL term) and "elevation property".

The Regulations ontology is built on top of a simplified version of ifcOWL, extended with complex concepts as pertaining to regulations. Regulation concepts are organized in layers. The ifcOWL vocabulary is the ground layer of the resulting vocabulary pyramid. At the very top of the pyramid, we find the regulation-specific vocabulary. A term of a specific layer is defined by using terms belonging to the layers beneath. This layered approach is flexible as it keeps high-level definitions simple by limiting the number of terms and using cascading rules among layers.

The process we implemented for providing a semantic regulatory view based on an IFC source model comprises several steps. First, we convert the IFC file into a RDF file using the Open-IFC-to-RDF-converter [START_REF] Pauwels | IFC-to-RDF Converter[END_REF] along with custom filtering of classes and relations. The output of this first step is stored in a repository of a triple store. Second, we use a "geometrical pre-processor" that renders and filters geometry aspects. Third, a "semantic preprocessor" filters and enriches the triple store knowledge base with terms from the Regulations ontology.

During this process, some IFC patterns are simplified [START_REF] Pauwels | SimpleBIM : from full ifcOWL graphs to simplified building graphs, Ework and ebusiness in architecture, engineering and construction[END_REF]. For instance, room classification is expressed in ifcOWL with at least ten triples, whereas the Regulations ontology allows narrowing this number down to one triple (see Fig. 1). Other IFC patterns are converted to higher level concepts from the Regulations ontology such as the "highest storey" concept.

Fig 1. Semantic shortcuts for simplifying classification assignment

From a regulatory point of view, the geometric pre-processor only passes significant classes and relations as highlighted in [START_REF] Preidel | Pre-Processing IFC Building Models for Code Compliance Checking based on Visual Programming[END_REF] and [START_REF] Daum | Processing of Topological BIM Queries using Boundary Representation Based Methods[END_REF]. Indeed, most of the time a bounding box is sufficient to express regulatory geometrical constrains, instead of relying IFC-based descriptions which contains large sets of tiny facets. Handling such constraints trough a bounding box simplifies their management in the system. Coordinates of each bounding box are expressed by triples.

From regulatory texts to query based constraints

This section provides details about the methodology used to transform the regulatory constraints, expressed in natural language, into processable queries.

First, regulatory texts are prepared and interpreted by working groups of experts. They rely on text editors' auto completion and syntax highlighting tools for re-writing regulatory texts into semi-formal rules. Each semi-formal rule is supposed to report building elements that are non-compliant with regard to a matter or regulation. Depending on its complexity, a regulatory rule expressed in a legal jargon can be divided into several atomic semi-formal rules dealing with complementary aspects. The idea is to keep each semi-formal rule as simple as possible. As a guideline, we suggested that each semi-formal rule begins with an "IF" statement followed by a condition on specific elements and ends with a "THEN NON-COMPLIANT" statement.

The regulatory constraints are pushed from semi-formal to formal stage. This step is performed by computer scientists familiar with BIM and Semantic Web technologies. When semi-formally expressed, rules are easier to adapt into formal rules. Each semiformal rule is translated into a SPARQL query using the same level of vocabulary coming from the Regulations ontology. Relying on such controlled vocabulary (e.g. the Regulations ontology) when expressing queries allows construction experts to still understand the underlying constraints and rules. Geometrical constraints are described by using geoSPARQL (OGC standard) concepts covering all topological relations that are possible between two geometries as introduced by McGlinn et al [START_REF] Mcglinn | Integrating Ireland's Geospatial Information to Provide Authoritative Building Information Models[END_REF]. Simplifying the geometry eases the handling of geometry coordinates. This approach should be extended if needed by using functional extensions, as it is the case in [START_REF] Zhang | Querying Linked Building Data Using SPARQL with Functional Extensions[END_REF].

Automating checking

Automating code-compliance checking consists of linking the output of conversion algorithm with the input of geometrical and semantic pre-processors as presented in section 3.1. Each time a new IFC model is submitted to the semantic checker this chain (algorithm + pre-processors) is triggered. The final output is a triple store populated with data with the correct level of details and aligned with the Regulations ontology. The geometrical pre-processor is java-based and returns geometrical triples from IFC. The semantic pre-processor consists of a sequence of semantic forward chaining operations that filters unnecessary information or enhances the knowledge base with high level vocabulary. Indeed, some IFC terms are almost never addressed by regulation texts -such as high-level geometry, element history or sensor states. Thus they can be removed from the knowledge base in order to reduce the amount of data. Ontology alignment and extensions to ifcOWL are declared in the Regulations ontology. SPARQL queries are organized by regulation topics (fire safety, accessibility, airing, acoustics…) so that the architect can focus on a specific aspect. Technically, constraint-queries are packed as SPARQL files in a ZIP archive. The archive also contains metadata for end-users' information purpose. Next steps for further automating algorithm consist in executing, one by one, each formal constraint-query corresponding to each chosen topic and verify the results.

The system outputs a Building Collaboration Format (BCF) file reporting for each constraint-query a list of non-compliant building elements. The BCF file makes it convenient to display such results within a standard IFC viewer. This allows displaying, on top of a graphical representation of a building element, the list of constraints that it breaks.

Results

We applied this approach on several constraints of the French regulation on both fire safety and accessibility domain. A first version of a French regulation ontology based on previous CSTB research works has been adapted [START_REF] Bouzidi | Towards a Semantic-based Approach for Modeling Regulatory Documents in Building Industry[END_REF]. A group of domain experts had analysed, interpreted and finally converted a dozen of regulation texts into about one hundred semi-formal constraints implementing concepts of the regulatory ontology vocabulary. The regulation ontology was then extended with experts' new concepts suggestions.

The following paragraphs illustrate in detail two specific use cases. According to French building regulation, a WC seat (water closet seat) must have a free space of 0.8 x 1 m at its left side OR at its right side. The "WC" is a high level regulatory concept corresponding to an IfcFlowTerminal with an IfcFlowTerminalType with a predefined type equal to "WCSEAT" [START_REF]IFC 2x3, Specification[END_REF]. The "FreeSpace" is also a high level regulatory concept. It represents a virtual object (not physical) with 3D dimensions (bounding box). This object shall not intersect any other physical element of the building. The checking process of this constraint (see Fig. 2) is detailed by the following sequence:

1. For each "WC" (water closet seat), two "FreeSpace" -and their bounding box-are created on-the-fly. One on the left and one on the right. 2. For each FreeSpace, the query retrieves all building elements that intersects. 3. Finally, the query returns all "WC" with at least one "FreeSpace" intersecting at least one building element.

In this use case, a bounding box representation of a building element is accurate enough to detect clashes with the "WC". Yet, the bounding box must be oriented to be able to locate left and right sides.

In addition to the list of elements that break rules, the explanation query displays the list of elements that intersects with a "FreeSpace". During our test we manage to detect a handrail that was not located on the right side of the "WC".

The fire safety domain regulation put some constraints on the structure element performances. Indeed, from the fire safety perspective, the structure elements of the building must remain bearing load during at least one hour in case of fire. This constraint differs according to the building height. Furthermore, the building height is measured from the ground until the floor of the latest storey. Translated as a SPARQL query, this constraint implies various high-level concepts such as: highest level (storey with highest elevation), storey floor (lower slab of a storey), structure element (building element that bear load), and loadbearing duration in case of fire. The following sequence is triggered by the constraint-query:

1. The building height according to the fire safety standard is calculated on the fly 2. The material performance threshold is determined according to the building height 3. All structure elements are retrieved and their performances are checked 4. Finally, the query returns all structure elements with fire resistance below the threshold.

This architecture has been implemented on the "French national BIM platform" providing a full automated checking process on various fire safety and accessibility topics.

Conclusion

As compared to MVDXML, we gained several advantages from our Semantic checker approach. At a first look, it needs transformation of an IFC model to the RDF representation which is an overhead and time taking process. But, once the semantic repository is prepared, it provides ease to unify, explore, analyze, and extend triples for the verification of IFC models. We conclude that semantic technologies provide more rich mechanisms and answer vast types of queries for the verification of IFC models. It provides powerful features based on SPARQL libraries and serves best for the automated code compliance and verification of IFC models.

While translating regulatory text to formal constraint-queries, we introduced a semiformal stage that is readable by both construction experts and data scientists. This study also points out that aligning the model to high level vocabulary makes it easier to write and maintain constraint-queries.

Semantic techniques make it possible to translate IFC terms to high level through two different way, i.e., forward chaining or backward chaining. Backward chaining consists of ontology statements that align IFC concepts with regulatory concepts, whereas forward chaining consists in creating supplementary triples. While, backward chaining using ontology statements is theoretically the best approach, we found that forward chaining is a good alternative at an optimization stage to save memory and CPU resources. From the machine point of view, backward chaining is processed each time a semantic query is submitted whereas forward chaining is executed each time the data changes.

We found that constraints included in regulations texts mostly consist in retrieving elements and checking their properties. Some constraints concerning geometry do not require a high level of details. Most frequent topological relations expressed in constraints deal with the intersection and adjacency. We also found that some information required when considering regulatory aspects are missing in the IFC specification. Some properties specific to national regulation need to be added to the IFC schema. Fire escape plan is also out of the IFC scope but is required when checking doors opening direction.

Obviously, this approach does not address all constraints defined in the regulations. For instance, when considering a complex corridor with multiples aisles -and their different related widths -the bounding box approximation will output erroneous results while computing intersection between spaces. Nevertheless, we argue that this approach can help in addressing the issues mentioned above. Furthermore, we envision extending it by considering the latest works done in the Linked Building Data Community group such as the Building Topology Ontology (BOT [START_REF] Rasmussen | Building Topology Ontology (BOT)[END_REF]), the ontology defining the core concepts of a building and OntoBREP [START_REF] Perzylo | An ontology for CAD data and geometric constraints as a link between product models and semantic robot task descriptions[END_REF] -ontology for CAD Data and Geometric Constraints.

Fig 2 .

 2 Fig 2. "FreeSpace around WC seat" constraint-query as a graph

Proceedings of the 6th Linked Data in Architecture and Construction Workshop