Nachum Dershowitz
email: nachum.dershowitz@cs.tau.ac.il

Jean-Pierre Jouannaud
email: jeanpierre.jouannaud@gmail.com

Drags: A Simple Algebraic Framework For Graph Rewriting

We are interested here in a natural generalization of term rewriting techniques to what we call drags, viz. finite, directed, ordered multigraphs. Each vertex is labeled by a function symbol, the arity of which governs the number of vertices to which it relates by an incoming edge in the graph. To this end, we develop a rich algebra of drags which generalizes the familiar term algebra and its associated rewriting capabilities. Viewing graphs as terms provides an initial building block for rewriting with such graphs, that should impact the many areas where computations take place on graphs.

Introduction

Rewriting with graphs has a long history in computer science, graphs being used to represent data structures, but also program structures, and even computational models. Since graphs originated from topological investigations, it is therefore no wonder to see them also play an important role in modern algebraic topology, more specifically in the theory of operads [START_REF] Dotsenko | Algebraic Operads: An Algorithmic Companion[END_REF]. Our interest in graphs originated from the latter theory which has strong relationships with various areas of computer science, including concurency theory and type theory.

Rewriting graphs is very similar to rewriting terms of an algebra, the same questions occur: What rewriting relation do we use? Is there an efficient pattern matching algorithm? Is rewriting terminating? Is it confluent?

The above questions have therefore been adressed by the rewriting community since the mid-eigthies. Termination and confluence techniques have been elaborated for various generalization of trees, such as rational trees, directed acyclic graphs, lambda-terms and lambda-graphs. But the design of a convenient mathematical framework for rewriting arbitrary graphs has made little progress beyond various categorical definitions and the study of the particular cases we just mentionned.

This paper describes the design of a general class of graphs -actually, multigraphs-that has good structural properties with respect to graph rewriting. In particular, graphs in this class can be naturally equipped with a natural, simple algebraic structure that allows to see them as terms with sharing, and therefore to mimic many techniques developped for trees in the past.

Despite being purely theoretical, the framework developped here should lead to easy implementations: the view of a drag as a term with sharing is extremely easy to implement, and therefore to convert many existing term rewriting implementations to actually do graph rewriting.

The particular class of graphs we are interested in is introduced in Section 2. The algebra of drag is described in Section 3, and the view of drags as terms in Section 4. Drag rewriting is introduced in Section 5.

Drags 2 Finite, directed, labeled multigraphs

The class of graphs with which we deal here is that consisting of finite directed multi-rooted graphs with labeled vertices and allowing multiple edges between vertices. In the present work, we assume that the outgoing neighbors (vertices at the other end of outgoing edges) are ordered (from left to right, say) and that their number is fixed, depending solely on the label of the vertex: we presuppose a set of function symbols Σ, whose elements f ∈ Σ used as labels are equipped with a fixed arity (we have no associativecommutative symbols here), and a denumerable set of variable symbols Ξ disjoint from Σ, whose arity is 0.

We shall call these finite directed rooted labelled graphs, drags.

The successors of a vertex labeled f in a drag are implicitly interpreted as the arguments of the function symbol f . The graph describes therefore the inverse of the computational flow. Some might prefer the converse choice. This class of graphs appears to be a very general model to describe computations in the absence of binding constructs.

To lessen notational burden, we use the vertical bars | • | to denote various quantities, such as length of lists, size of sets or of expressions, and even the arity of function symbols. We use ∅ for an empty list, set or multiset, ∪ for list concatenation as well as set and multiset union, \ for list, set or multiset difference, and identify a singleton list, set or multiset with its content.

In contrast with the standard categorical approach [START_REF] Ehrig | Graph-Grammars: An Algebraic Approach[END_REF][START_REF]Handbook of Graph Grammars and Computing by Graph Transformation[END_REF], our multi-graph model is very standard. The novelty, however, is that we allow for arbitrarily many roots and and arbitrarily complex cycles. Definition 1 (Drag). A (closed) drag is a tuple V, R, L, X , where V is the finite set of vertices; R is a finite list of vertices, called roots, such that R(n) is the nth root in the list; L : V → Σ is the labeling function, mapping vertices to labels from some vocabulary Σ; and X : V → V * is the successor function, mapping each vertex v ∈ V to a list of vertices in V whose length equals the appropriate arity |L(v)|. Given b ∈ X(a), then (a, b) is an edge with source a and target b. We also write aXb. The reflexive transitive closure X * of the relation X is called accessibility. A vertex v is said to be accessible if rX * v for some root r. A drag is clean if all its vertices are accessible. A root r is maximal if ∀r ∈ R. r X * r implies rX * r .

Labelings extend to lists, sets and multisets of vertices as expected.

Definition 2 (Open drag

). An open drag is a drag over the set Σ ∪ Ξ, where Ξ is a set of variables. The vertices labeled by a function symbol in Σ are internal; those labeled by a variable are called sprouts.

When non-empty, the set S of sprouts will be added at the end of the tuple which becomes V, R, L, X, S . An open drag is linear if different sprouts have different labels. It is cyclic if R∩S = ∅ and ∀v ∈V \S. ∃r ∈ R. vX * r.

The component R is a list with repetitions, whereas S is a set. Allowing for repeated roots is needed for rewriting. Sprouts may be roots. Cleanness relates to garbage collection, discussed further in Section 5. It will sometimes be convenient to consider roots as specific incoming edges and to identify a sprout with its label.

A cyclic drag may consist of a single internal vertex being a root. Although this may seem odd to call such a drag cyclic, this is an essential ingredient here which allows to see a tree as a particular drag. For example, the ground term a and the term f (x) are cyclic drags, while the terms f (a) and f ( f (x)) are not. Terms, sequences of terms and terms with sharing are particular drags. Nonlinear drags play an essential rôle for sharing and unsharing subterms in terms.

Given a drag D, we use the following notations: V er(D) for its set of vertices; X D for its successor function; A cc(D) for its set of accessible vertices; R(D) for its set of roots; R(D) max for the subset of maximal roots; S (D) for its set of sprouts; V ar(D) for the set of variables labeling its sprouts; D r s to indicate its roots and sprouts; and D ↓ for the clean drag obtained by removing its inaccessible vertices and related edges. Given finally L,V which are respectively a list and a set of vertices of some drag D, we denote by L \V , the sublist of L obtained by filtering out the vertices belonging to V .

Subdrags

Definition 3 (Subdrag). Given a drag U = V, R, L, X, S , and a subset W of its vertices, the subdrag U |W of U generated by W is the drag V , R , L , X , S where (i) V is the least superset of W that is closed under X; (ii) L , X , S are the restrictions of L, X, S to V ;

(iii) R is (R ∩V ) ∪ (X(V \V ) ∩V ).
A subdrag is clean by construction. Its roots are obtained by adding as new roots those vertices which have an incoming edge in U that is not in U |W . The order of elements in this additional list does not really matter for now.

In particular, restricting a drag D to its maximal roots, yields D ↓. On the other hand, restricting any drag to an empty set of vertices, yields the empty drag

∅ def =1 ∅ ∅ = (∅, ∅, ∅, ∅, ∅).

Drag Isomorphisms

Drags are graphs. An isomorphism between two open drags is a one-to-one mapping between their respective sets of (accessible) vertices that identifies their respective labels (up to renaming of the sprouts' labels done here by the very same mapping) and lists of roots, and commutes with their respective successor functions. It is sometimes useful to consider multisets instead of lists of roots, resulting in a coarser equivalence on drags:

Definition 4 (Drag isomorphism). Given two open drags D = V, R, L, X, S and D = V , R , L , X , S whose sprouts are identified with their labels, an isomorphism from D to D is a one-to-one mapping o : A cc(V ) → A cc(V ) such that (i) o restricts to bijections between the lists of roots and sets of accessible sprouts;

(ii) ∀v ∈ A cc(V ) \ S. L(v) = L (o(v)) and X(v) = X (o(v));
We Note that inaccessible vertices are ignored, making any drag D isomorphic (with ≡) to D ↓. Identifying the sprouts with the variables that label them saves an additional one-to-one mapping between the variables that would otherwise become necessary.

Example 5. Our running example is made of seven single-rooted drags represented in Figure 1.

The pictorial representation does not tell us precisely which argument of a function symbol comes first. In our convention, the leftmost one in the figure comes first. Another convention tells us the order of arguments, by sweeping the plane counterclockwise from the first. Finally, an incoming down-arrow indicates a root.

Ariola and Klop's wording [START_REF] Ariola | Cyclic Lambda Graph Rewriting[END_REF], horizontal/vertical sharing, is intuitive.

Let Σ = {a, b, f , g, h} with | f | = 3, |h| = 2, |g| = 1, |a| = |b| = 0. We first consider the term f (g(b), a, g(b)) in which b and g(b) can Drags ↓ f © g c b c a d d g c b ↓ f © g c a d d g d d © b ↓ f © g c b c a ↓ f © a c g c b ↓ f © g c g d d h d d © b ↓ f © g g c h © h c c g g ↓ f © b← c h d d h c g e e u
Figure 1: Seven non-isomorphic drags whose vertices are identified with the extremities of their edges. be shared. Among all possible cases, three are represented. The 5th has both horizontal and vertical sharing. The 6th has three cycles: an independent leftmost loop and a horizontal cycle sharing another horizontal cycle below. The last has a single big cycle, which is both horizontal and vertical including an inside cycle.

Drag algebra

We equip drags with a composition operator, which, given two drags, connects the sprouts of each drag with the roots of the other. As usual, we denote by Dom(ξ ) and Im(ξ ) the domain (of definition) and image of a (partial) function ξ , using ξ A→B for its restriction going from subset A of its domain to subset B of its image, omitting →B when irrelevant. Both are sets. hence Im(ξ ) may differ from ξ (R), which is a list with possibly repetitions when R is a list itself. First, comes the connection device: Conditions (i,ii) are always satisfied for linear drags, or for nonlinear drags whose switchboard, called linear, is defined for sprouts whose variables are all different. When this is not the case, checking that a particular ξ is a switchboard involves graph isomorphism, which is polynomial in our case, since we have ordained fixed arities. (Were one to allow variadic vertices, isomorphism would be quasipolynomial [START_REF] Babai | Graph Isomorphism in Quasipolynomial Time[END_REF].) A simple effective approximation here is to require instead that R (ξ D (s)) = R (ξ D (t)) for (i) and similarly R(ξ D (s)) = R(ξ D (t)) for (ii). This approximation is not expressive enough for the case of trees, but turns out to fit shared trees, for which isomorphism is computable in linear time. In the sequel, we assume this latter condition for simplicity.

Definition 6 (Switchboard). Let D = V, R, L, X, S and D = V , R , L , X , S be open drags. A switch- board ξ for (D, D ) is a pair ξ D : S → [1 .. |R |], ξ D : S → [1 .. |R|]
In practice, we shall usually enforce V ar(D) ∩ V ar(D ) = ∅ so as to avoid confusion, but this is by no means necessary since switchboards are pairs of mappings operating on V ar(D) and V ar(D ) independently, not a single mapping operating on V ar(D) ∪ V ar(D ). Further, splitting ξ makes the definition clearly symmetric.

We shall often restrict switchboards to subdrags, denoting by ξ A→B , the restriction of ξ D going from the subdrag A of D to the subdrag B of D , using ξ A when B = D .

Note finally that injectivity of ξ and the fact that roots are lists with repetitions go along together. One could think of using sets for roots, in which case ξ would not be injective. This alternative, however, does not fit well with rewriting, for which it will be essential to control precisely both the number and order of roots of left-and right-hand sides of rules. Further, injectivity implies that the list ξ D Dom(ξ D ) is indeed a set, hence making the set difference

(V \ Dom(ξ D )) ∩ (V \ Dom(ξ D )) = ∅. Their composition is the drag D ⊗ ξ D def = V , R , L , X , S , where (1) V = (V \ Dom(ξ D )) ∪ (V \ Dom(ξ D )); (2) R = R([1..|R|] \ ξ D (S )) \ Dom(ξ D ) ∪ R ([1..|R |] \ ξ D (S)) \ Dom(ξ D ) ∪ R(ξ D (S ∩ (R \ R (ξ D (S))))) ∪ R (ξ D (S ∩ (R \ R(ξ D (S ))))); (3) L (v ∈ V \Dom(ξ D )) = L(v) and L (v ∈ V \Dom(ξ D )) = L (v); (4a) X (v ∈ V \ S) = if X(v) ∈ Dom(ξ D ) then X(v); otherwise R(ξ D (v)); (4b) X (v ∈ V \ S ) = if X (v ) ∈ Dom(ξ D ) then X (v ); otherwise R (ξ D (v )); (5) S = (S \ Dom(ξ D )) ∪ (S \ Dom(ξ D )).
The essence of this definition is that the (disjoint) union of the two drags is formed, but with sprouts in the domain of the switchboards merged with the vertices referred to in the switchboard images. This necessitates renaming (via ξ ) the targets of edges that had pointed to the sprouts. The only complication is the calculation of the list of roots R of the resulting drag D . This list is made of two sublists for each component drag, we describe them for D:

1. the remaining roots of D: they are obtained by removing from R the roots which are in the image of ξ D (note that we remove the root indices from their set before to apply R), before to eliminate all those sprouts belonging to the domain of ξ D which are also roots.

2. the vertices of D belonging to the image by ξ D of a sprout of D which was a root in D : vertices inherit the root status from the sprouts with which they are merged. Of course, one needs to worry about the case where multiple sprouts are merged successively, when the switchboards map sprout to rooted-sprout to rooted-sprout, possibly in a cycle. In such cases, all need to be renamed to the same one, while accumulating the root statuses. Note that if ξ is a rewriting switchboard, then all roots and sprouts of D disappear in the composed drag. Otherwise, the symmetry of the definition is broken by choosing the roots originating from D to come first.

It is easy to see that D ⊗ ξ D is a well-defined drag, that is, it satisfies the arity constraint required at each vertex. Example 9. We show in Figure 2 three examples of compositions, the first two with similar drags. The first composition uses a directed switchboard, while the second uses a rewriting switchboard which induces a cycle. Variables x, y label the sprouts. In the second example, the remaining root is the first (red) root of the first drag which has two roots, the first red, the other black. The third shows how sprouts which are also roots connect to roots in the composition. Our definition of switchboard being (almost) symmetric, composition is itself symmetric provided all roots of the resulting drag originate from the same side. Otherwise, composition yields drags that are equal up to a cyclic permutation of their roots.

Drags ↓ f ↓ x ⊗ {x →1} ↓ f ↓ y = ↓ f ↓ f ↓ y ↓↓ f ↓ x ⊗ {x → 1, y → 2} ↓ f ↓ y = ↓ f ↓↑ f 1 ↓ f 2 ↓ f x ⊗ {x →1} 1 2 ↓ a = 1 ↓ f 2 ↓ f 3 ↓ a
Lemma 10 (Commutativity). Composition of drags is quasi-commutative, and becomes commutative for switchboards whose context or substitution is surjective, e.g. rewriting switchboards.

Commutativity shows that there is no big difference between the context and the substitution of a switchboard in a composition D ⊗ ξ D , since they simply exchange each other when exchanging the input drags of the composition.

Let us call identity a linear open drag whose all vertices are its sprouts and set of edges is empty. We denote it by 1 Y Z , where Y ⊆ Z * is its list of roots. We use ∅ for the drag 1 ∅ Z , which is isomorphic to 1 ∅ ∅ , and is called the empty drag for that reason.

Lemma 11 (Neutral). Let D be an open drag, X ⊆ V ar(D), and ι the directed identity switchboard for (D, 1 X X ). Then D ⊗ ι 1 X X = D. In particular, this property holds when X = ∅.

Definition 12 (Compatibility). Two switchboards ξ and ζ for the respective pairs of drags (U,V ) and

(V,W ) are compatible if (i) Dom(ξ V ) ∩ Dom(ζ V ) = ∅ and (ii) Im(ξ U ) ∩ Im(ζ W ) = ∅.
Lemma 13 (Associativity). Let U,V,W be three drags, and ξ , ζ be compatible switchboards for (U,V ) and (V,W ) respectively. Then,

(U ⊗ ξ V ) ⊗ ζ W = U ⊗ ξ (V ⊗ ζ W ).
Proof. Compatibility ensures that ξ and ζ are switchboards for (U,V ⊗ ζ W ) and (U ⊗ ξ V,W ), by eliminating the possibility that a sprout or root is used twice. Both sides of the associativity equation are thus defined. The equality itself results from easy checking.

Drags decomposition

Next, we investigate ways of decomposing a drag. We first give a general construction which splits a drag into the subdrag generated by some of its vertices, and the rest of the drag, its antecedent. Lemma 14. Given a drag D and a subset of vertices W ⊆ V , there exists a linear drag A, called antecedent, and a total, surjective, and directed switchboard

ζ such that D = A ⊗ ζ D |W .
The fact that the switchboard ξ is directed expresses the property that decomposing a drag into a subdrag and its antecedent does not break any of its cycles. Note that ξ induces an order on the roots of the subdrag which are not roots of the whole drag. If W =V , then (A, ξ ) must be the identity extension 1 X X with |X| = |R(D)|.

A tree headed by the symbol f of arity n has one root labeled by f , or equivalently a head corresponding to the expression f (x 1 , . . . , x n ), and several subtrees, seen here as a single drag with several

1 2 3 ↓ g c b ↓ a ↓ g c b 1 2 3 ↓ g ↓ a ↓ g d d © b 1 3 2 ↓ g c b ↓ a 1 2 3 ↓ g c b ↓ a 1 2 3 ↓ g ↓ g ↓ h d d © b 1 2 3 ↓ g ↓ h h c c g g 1 2 3 ↓ b← ↓ h ↓ h c g e e u 1 2 ↓ b ↓ b 1 2 ↓ b 1 ↓ b 1 ↓ b 1 2 ↓ b 1 2 ↓ g ↓ g 1 2

↓ b

Figure 3: Successive subdrags of the drags of Figure 1 roots. This unique decomposition is a fundamental property of trees which expresses the fact that the set of trees equipped with the "head-operations" has an initial algebra structure. Likewise, a drag has a head, which is intended to be the largest cycle containing the maximal roots of the drag -also the smallest nontrivial antecedent of the drag -and one tail, possibly a list of several connected components. The cycle may be of zero length in order to have trees as a particular case. As for trees, the head will have a list of new sprouts which are in one-to-one correspondence with the roots of the subdrag. A drag can therefore be seen as a bipartite graph made of its head, its tail, and a directed switchboard specifying that correspondence.

Definition 15. The tail ∇D of a drag D = V, R, L, X, S is the subdrag generated by the set of vertices

V \ {v ∈ V : vX * R(D) max }.
As an application of Lemma 14, we get: Lemma 16. Given a drag D, there exists a linear drag D, called head of D, and a directed switchboard

ζ such that D = D ⊗ ζ ∇D.
The head of a drag is therefore the antecedent of its tail. Note that ξ is a bijection between the sprouts of the head and the roots of the tail that were not already roots of the drag D. In the sequel, we assume that the sprouts of the head are ordered with respect to a depth first search initialized with its list of roots. This order on the sprouts of the head induces via the switchboard ξ , an order on the roots of the associated tail which were not roots of the original drag. The tail is therefore now canonically determined.

We could have defined the head and tail of D as a pair of drags whose head contains all maximal roots, is cyclic, and its composition with the tail by a directed switchboard is the drag D itself. This characterization of the head of a drag is a property that we should keep in mind.

We can now show the structure theorem of open drags:

Theorem 17 (Structure). Isomorphic drags have isomorphic heads and tails.

Proof. Given two isomorphic drags D, D , constructing bijections between their respective heads and tails from the bijection between their sets of vertices is an easy exercise.

Example 18. The successive tails of the drags of Figure 1 are represented in Figure 3. The numbering displayed above the roots of each tail specifies the order of roots in their list of roots. These tails, represented on the first row, have all the same number of roots which is determined by their head f . The next row represents the next level of tails, which here have one or two roots, depending on their own head. The last tail in that row has two roots, a nonmaximal one in its father drag, and a new one. Identity extensions allow one to change the structure of a drag without changing its internal nodes. If the drag has a single root, it is easy to see that identity extensions are enough to predict all forms that a drag may take under composition with an extension. This is no longer true with multi-rooted drags, since identity cyclic extensions cannot reach two different roots from the same sprout. Note that non-identity cyclic extensions cannot do any better if the vocabulary contains unary symbols only, as in Example 20. If the vocabulary contains a symbol g with |g| = 2, then new cyclic extensions pop up, both with the same drag ({1, 2, 3}, {1}, L(1) = g and L(2) = y and L(3) = z, X(1) = {2, 3}, {2, 3}), and with the respective

1 ↓ f 2 3 ↓ f x L = L ⊗ ∅ 1 ∅ ∅ 2 3 ↓ f f L ⊗ {x →y,y →1} 1 {y} {y} 1 ↓ f 2 f L ⊗ {x →y,y →2} 1 {y} {y}
switchboards {x → 1, y → 1, z → 2} and {x → 1, y → 1, z → 3}.
Here is the property at the heart of drag rewriting:

Lemma 21 (Decomposition). Let U,W be clean nonempty drags and ξ be a switchboard for (W,U). If ξ W is surjective on R(U) max , then, there exist drags A, B,C and switchboards ζ , θ such that 1 (B, ξ B , ξ U→B ) is a cyclic extension of U denoted by (B, ξ );

2 (C, θ ) is a substitution extension of B ⊗ ξ U; 3 (A, ζ ) is a context extension of (B ⊗ θ U) ⊗ η C; 4 W ⊗ ξ U = A ⊗ ζ ((B ⊗ θ U) ⊗ η C);
5 C is empty if all internal nodes of W reach one of its sprouts.

Proof. The lemma is depicted in This lemma applies of course to rewriting extensions (W, ξ ) of U, since ξ W is surjective on all roots. It applies as well to cyclic extensions (B, ξ ) since ξ U ) generates B, hence must be surjective on its maximal roots. but Definition 22 (Rules). A graph rewrite rule is a pair of open clean drags written L → R such that |R(L)| = |R(R)| and V ar(R) ⊆ V ar(L). A graph rewrite system is a set of graph rewrite rules.

Trees are particular clean drags. Because they have a single root, term rewrite rules satisfy the first condition. The second must be explicitly stated in the definition of a term rewrite rule, as for drags.

Note that we could allow R to be unclean, and strengthen instead the second condition to be V ar(R) = V ar(L). Adding unreachable sprouts to a clean R would then do with an isomorphic right-hand side. Such a handy trick is not possible with trees. All assumptions on ξ play a rôle. First, because ξ is a rewriting switchboard, ξ W must be linear. This implies that the variables labeling the sprouts of W which are not already sprouts of D must be all different. Second, ξ C must be surjective, implying that the roots of L disappear in the composition. Third, ξ L must be total, implying that the sprouts of L disappear in the composition. Fourth, D must be nonempty, implying that the roots of W do not all disappear in the composition; hence ξ L cannot be surjective. Finally, if L is linear, then ξ is a linear rewriting switchboard: no test for isomorphism is needed when pattern matching D with L.

Drags

A B U C ζ ξ A ζ ξ B ξ U→C θ v t r s s v ξ U→B
If left-hand and right-hand sides of a drag rewrite rule have the same variables as sprouts, the computed drag generated will be clean if the starting drag is. But if there are strictly fewer variables in the right-hand side, replacement may generate inaccessible vertices, hence require garbage collection, which is therefore explicitly built in this model of computation. On the other hand, duplication is banned by the drag rewriting model, which therefore generalizes shared rewriting, but not term rewriting per se.

Can our drag model capture term rewriting ? As should be expected, the anwer is affirmative. The term model does not permit sharing. Taking terms with sharing as our term structure, terms without sharing become normal forms in that richer structure. We can therefore capture term rewriting within our model by changing the definition of a clean drag. Details are left to the reader.

In the case of trees, the rewrite extension of L is made of a context and a substitution. What is remarkable about the above definition is that there is no longer any distinction between the two, which, together, constitute the rewriting extension of the left-hand side. The context and substitution of L correspond respectively to the context extension and substitution extension of Lemma 21. The cyclic extension is of course the identity in the case of trees.

Example 24. Consider the cyclic composition of Figure 2, and let f (y ) → y be a graph rewrite rule whose left-hand and right-hand sides are, more precisely, the respective drags

({1, 2}, 1, L(1) = f and L(2) = y , X(1) = 2, {2}) and ({1}, 1, L(1) = y , , {1}). Then, the drag ({1, 2}, 1, L(1) = L(2) = f , X(1) = 2 and X(2) = 1) rewrites to the drag ({1, 2}, 1, L(1) = f and L(2) = x , X(1) = 2, {2}) ⊗ {x →1,y →1} ({1}, 1, L(1) = y , , {1}) = ({1}, 1, L(1) = f , X(1) = 1): a cycle of length 2 has been rewrit- ten into one of length 1.
It is also possible to break a cycle in a drag: the same drag ({1, 2}, 1,

L(1) = L(2) = f , X(1) = 2 and X(2) = 1) rewrites to the drag ({1, 2}, 1, L(1) = f and L(2) = x , X(1) = x , {2}) ⊗ {x →1} ({1}, 1, L(1) = a, ) = ({1, 2}, 1, L(1) = f , L(2) = a, X (
↓ f ↓↑ f = 1 2 ↓ f ↓ x ⊗ {x → 1, y → 2} 1 ↓ f ↓ y -→ 1 2 ↓ f ↓ x ⊗ {x → 1, y → 2} 1 ↓ y = 1 ↓ f (shrinking) 1 ↓ f ↓↑ f = 1 2 ↓ f ↓ x ⊗ {x → 1, y → 2} 1 ↓ f ↓ y -→ 1 2 ↓ f ↓ x ⊗ {x → 1, y → 2} 1 ↓ a = 1 ↓ f ↓ a (breaking)
Figure 6: Rewriting a cycle.

These rewrites are shown in Figure 6. In both cases, the upper occurrence of f (in blue) is part of the context, while the one below (in red) is part of the rewrite rule. Rewriting the upper occurrence of f instead of the lower is an exercise left to the reader.

Lemma 25. If U -→ R V , then V ar(V ) ⊆ V ar(U).

Proof. The property is true of rules and preserved by rewriting.

Lemma 26. Assume U,V,W are three open drags such that the pair (U,V ) is a rewrite and ξ is a switchboard for (W,U). Then, ξ is a switchboard for (W,V ).

Proof. Because the lists of roots of U and V have the same length, and because the sprouts of U,V are labeled by the same variables that are here identified with the sprouts themselves, ξ U does need any change, and since V ar(V ) ⊆ V ar(U), ξ V is the restriction of ξ U to the sprouts of V whose variables are in V ar(W ).

A consequence of this property is that properties of extensions are preserved by rewriting. Further, the decomposition lemma itself is also preserved:

Conclusion

We have invented drags, a very general class of multigraphs, and explored their nice and useful properties, which yield a new, simple and very effective model for graph rewriting. In a companion paper, we have used this model to show that termination techniques that have been developped for terms, namely the recursive path order, scale to drags, but we believe that this is a general fact for syntactic techniqueds. This allowed us, in particular, to design a total order for extending the theory of Gröbner bases to algebraic operads, which are just polynomial expressions built over drags.

Our term representation and ordering have the potential to be further extended with abstraction and application, as done in the computability path order [START_REF] Blanqui | The Computability Path Ordering[END_REF], and obtain an ordering over λ -terms with sharing and back-arrows, also called lambda graphs [START_REF] Ariola | Cyclic Lambda Graph Rewriting[END_REF]. We presume that this problem conceals some unseen difficulties.

  write D = o D , or simply D = D when D and D are isomorphic drags. We write D o D , or simply D D , and say that D, D are quasi-isomorphic if o restricts to a bijection between the multisets (instead of lists) of roots. We write D ≡ D and D ∼ = D instead of D = D and D D in case o is the identity.

  of partial injective functions (we also say that (D , ξ ) is an extension of D), such that (i) ∀s,t ∈ S. s ∈ Dom(ξ D ) and L(s) = L(t) imply t ∈ Dom(ξ D ) and D |R (ξ D (s)) D |R (ξ D (t)) ; (ii) ∀s,t ∈ S . s ∈ Dom(ξ D ) and L (s) = L (t) imply t ∈ Dom(ξ D ) and D |R(ξ D (s)) D |R(ξ D (t)) .

  [1..|R |] \ ξ D (Dom(ξ D )) well defined. Two particular kinds of switchboards play an important rôle: Definition 7 (Categorization). A switchboard ξ for (D, D ) is • directed if one of ξ D and ξ D has an empty domain; the pair (D , ξ D ) is called a context extension of D if Dom(ξ D ) = ∅, and a substitution extension of D if Dom(ξ D ) = ∅; • a rewriting switchboard if ξ D is linear and surjective and ξ D is total; the pair (D , ξ ) is called a rewriting extension of D. Directed switchboards correspond to the tree case, with all connections from one drag to the other. Rewriting switchboards allow one to "encompass" drag D in drag D. Note that context extensions and substitution extensions exchange each other by symmetry of the definition of a switchboard. A switchboard induces a binary operation on open drags: Definition 8 (Composition). Let D = V, R, L, X, S and D = V , R , L , X , S be open drags, and ξ be a switchboard for (D, D ) such that

Figure 2 :

 2 Figure 2: Unidirectional, cyclic and root-transfer compositions.

Figure 4 :

 4 Figure 4: Minimal cyclic extensions.

Figure 5 .

 5 Let first D be the subdrag of W ⊗ ξ U generated by R(U) max . Since U is clean, D must contain all vertices of U by definition of R(U) max . By Lemma 14, W ⊗ ξ U = A ⊗ ζ D for some context extension (A, ζ ) of D such that ζ is total, surjective and directed. Observe that ζ A→U = ξ A . Let now C be the subdrag of D generated by those vertices which cannot reach a sprout in Dom(ξ W ), hence a root of U. By definition of C and the fact that D contaions all vertices of U, the ancestor of C in D must contain all vertices of U, hence is of the form B ⊗ ξ B ,ξ U→B U. By Lemma 14, D = (B ⊗ ξ U) ⊗ θ C for some θ which is total, surjective and directed. Observe that θ U = ξ U→C . Property (4) now follows. Property (5) is straightforward.

Figure 5 :

 5 Figure 5: Decomposition of a pair of drags.

  1) = 2) with the rule f (y ) → a whose both sides are the drags ({1, 2}, 1, L(1) = f and L(2) = y , X(1) = 2, {2}) and ({1}, 1, L(1) = a, ).

1

 1

Lemma 27 .

 27 Assume U-→ R V . Then, U = A ⊗ ζ ((B ⊗ ξ L) ⊗ θ C) and V ≡ A ⊗ ζ ((B ⊗ ξ L) ⊗ θ C)for some A, B,C and ζ , ξ , θ such that (A, ζ ), (B, ξ ) and (C, θ ) are context, cyclic and substitution extensions, resp. Proof. Using Lemma 21.

Drags

The fundamental property of a drag expressed by Theorem 17 is that its decomposition into a head and tail is a faithful representation of it. This property holds true because drags are multi-rooted. Unirooted drags cannot represent horizontal sharing, in contrast to vertical sharing which can sometimes be preserved with uni-rooted drags. Moreover, the fact that tails are quasi-isomorphic does not hamper faithfulness: different orders of new roots for the tails yield different switchboards but the order of roots resulting from the composition does not depend upon those orders.

Drag rewriting

The notion of composition lead in the previous section to an initial, canonical way of decomposing a drag, into its head and tail. In this section, we investigate further ways of decomposing a drag, possibly breaking its cycles, which will of course become important when rewriting drags. More precisely, can a drag D be seen as the composition of a given drag U with some context W via some switchboard ξ ? In this case, we would of course say that D matches U, W and ξ being the matching context and switchboard.

The idea is that W splits into three parts: the vertices of W that are accessible from some sprout of U and can access some of its roots define a drag B; those that are accessible from some sprout of U but cannot access any of its roots define a drag C which is therefore a substitution extension of U; those which are not accessible from the sprouts of U form the remaining part A, which is therefore a context extension of U.

Particular substitution extensions of U are identity extensions

where x 1 , . . . , x n are variables labeling sprouts of U, and, of course, z ∈ Dom(ξ ). Then, the drag resulting from the composition of U with that extension is the same as U, except that all its sprouts labeled by the variables x 1 , . . . , x n are now merged into a single sprout labeled by z. The rôle of these identity extensions is to modify the structure of U by introducing sharing among its sprouts.

We are left defining the kind of extension that B belongs to: Definition 19 (Cyclic extension). An extension (B, ξ ) of a clean drag U is cyclic if B is a linear drag generated by Im(ξ U ), B⊗ ξ U is a clean nonempty drag, and for all t ∈V , there exists s ∈ Dom(ξ B ) such that tX * s. The extension is trivial if B is an identity drag and total if ξ B is surjective.

This condition resembles the one for cyclic drags, but does not imply that B ⊗ ξ U is cyclic because not all internal vertices of U may reach a root of B ⊗ ξ U. Only total extensions are cyclic, up to possibly ground subdrags of U when there are such subdrags in U.

The conditions for being a cyclic extension impose that ξ U is surjective on R\S as a set so as to generate B. The roots of the resulting drag are therefore either remaining roots of B (possibly some roots in S transfered to B ⊗ ξ U) or roots of U, but there must be sufficiently many of them so that the resulting drag is clean.

Identity cyclic extensions of U are of the form (1 Z Y , ι), where the variables in Y are one-to-one with those in Dom(ι U ) ⊆ V ar(U), and ι 1 : Y → [1 .. |R|] is an arbitrary map. The rôle of cyclic extensions is to modify the structure of U by connecting some of its sprouts to some of its roots. Unfortunately, identity cyclic extensions do not suffice for that purpose unless in special cases: Example 20. Let L be the be the drag made of two copies of the tree f (x) sharing the variable x. Three identity extensions of L are represented at Figure 4. The first is the trivial directed extension while the two others are cyclic ones. The second maps the variable x of L to the only root of the identity drag, which is its sprout y, and the only sprout y of the identity drag to the first root of L. The third maps instead y to the second root.


Pablo Alvarez
email: pablo.alvarez@etudiant.univ-rennes1.fr

Matthieu Chabanas

Simon Rouzé

Miguel Castro

Yohan Payan

Jean-Louis Dillenseger

Lung deformation between preoperative CT and intraoperative CBCT for Thoracoscopic Surgery: a case study

Keywords: Video Assisted Thoracoscopic Surgery (VATS), Image Registration, Lung, Cone-Beam CT

Video-Assisted Thoracoscopic Surgery (VATS) is a promising surgical treatment for early-stage lung cancer. With respect to standard thoracotomy, it is less invasive and provides better and faster patient recovery. However, a main issue is the accurate localization of small, subsolid nodules. While intraoperative Cone-Beam CT (CBCT) images can be acquired, they cannot be directly compared with preoperative CT images due to very large lung deformations occurring before and during surgery. This paper focuses on the quantification of deformations due to the change of positioning of the patient, from supine during CT acquisition to lateral decubitus in the operating room. A method is first introduced to segment the lung cavity in both CT and CBCT. The images are then registered in three steps: an initial alignment, followed by rigid registration and finally non-rigid registration, from which deformations are measured. Accuracy of the registration is quantified based on the Target Registration Error (TRE) between paired anatomical landmarks. Results of the registration process are on the order of 1.01 mm in median, with minimum and maximum errors 0.35 mm and 2.34 mm. Deformations on the parenchyma were mesured to be up to 14 mm and approximately 7 mm in average for the whole lung structure. While this study is only a first step towards image-guided therapy, it highlights the importance of accounting for lung deformation between preoperative and intraoperative images, which is crucial for the intraoperative nodule localization.

INTRODUCTION

Lung cancer remains as the worldwide leading cause of cancer death for both women and men. [START_REF] Jemal | Global Cancer Statistics[END_REF][START_REF]World Cancer Report 2014[END_REF] Such a high mortality is related to the late detection of the disease, where curative treatements are normally not available and the 5-year survival rate lies between 6% and 18%. [START_REF]World Cancer Report 2014[END_REF][START_REF] Siegel | Cancer Statistics[END_REF] However, screening programs performed on patients at risk have demonstrated that survival rates might be significantly increased if diagnosis and treatement are performed at early stages. [START_REF] Henschke | Early Lung Cancer Action Project: overall design and findings from baseline screening[END_REF][START_REF] Lung | Reduced Lung-Cancer Mortality with Low-Dose Computed Tomographic Screening[END_REF] In such scenarios, surgical resection of malignant nodules is prescribed to patients. The treatement is performed via either open thoracotomy or video-assisted thoracoscopic surgery (VATS), the latter being the least invasive method with better and faster patient recovery. [START_REF] Mckenna | New approaches to the minimally invasive treatment of lung cancer[END_REF] Even if preoperative CT images are used for planning VATS intervention, intraoperative localization of lung nodules is still challenging in many cases. This is particulary true when the nodules to be resected are small, sub-solid or deep within the parenchyma. [START_REF] Gould | Evaluation of Individuals With Pulmonary Nodules: When Is It Lung Cancer?[END_REF] The cause of this problem is the anatomical disparity between the intraoperative and preoperative configurations, as a consequence of large lung deformations present during surgery. These lung deformations can be mainly associated to two different sources: on the one hand, the patient position is changed from supine in preoperative CT acquisition to lateral decubitus in the operating room, which affects the way gravity influences internal organs. On the other hand, for the comfortable manipulation of the lung during surgery, the surgeon creates space inside the thoracic cage by allowing air getting into the intrapleural space. This phenomenon, known as a pneumothorax, produces a total collapse of the lung towards the mediastinum that modifies internal lung structures.

Thanks to its low dose radiation and fast acquisition time, intraoperative Cone-Beam CT (CBCT) imaging could guide the localization of challenging nodules during a VATS intervention. [START_REF] Uneri | Deformable registration of the inflated and deflated lung in cone-beam CTguided thoracic surgery: Initial investigation of a combined model-and image-driven approach[END_REF][START_REF] Rouzé | Small pulmonary nodule localization with cone beam computed tomography during video-assisted thoracic surgery: a feasibility study[END_REF] Nonetheless, lung structures are more difficult to see in CBCT images given the reduced image quality when compared to CT images. In addition, the intensity contrast between lung nodules and lung's parenchyma is decreased as a consequence of the deformation induced by the pneumothorax. In fact, the parenchyma becomes denser due to lung deflation. This is particularly problematic for the localization of low density lung nodules also referred as Ground Glass Opacities (GGO). A possible solution to this problem might be the superposition of preoperative CT information (e.g. segmentation of nodules and other important structures) with the intraoperative CBCT image, via an image registration procedure. However, the existence of large lung deformations makes such a task a real challenge. This paper focuses on the quantification of the deformations induced by the change of the patient position between preoperative and intraoperative configurations during a VATS intervention. Understanding these deformations might be an important factor towards the development of an efficient image-guided surgery procedure. A non-rigid registration method is proposed for the superposition of preoperative CT and intraoperative CBCT lung structures. The resulting geometrical transformation is then used to quantify the deformations needed to achieve such superposition. To the best of our knowledge, no study has ever tried to quantify the lung deformations occurring during a VATS intervention, even before pneumothorax.

MATERIALS AND METHODS

Data

The work herein presented is a feasibility study performed on one clinical case only. In this context, a wedge resection was prescribed to the patient for a solitary nodule of approximately 13 millimeters in diameter. The VATS surgical intervention was done at the Rennes University Hospital, Rennes, France. This study was approved by the local ethics committee and the patient gave informed consent prior to the procedure.

The study consisted of the acquisition of two tomographic images: a preoperative CT, where the patient is in supine position and was instructed to hold his breath during the capture (end of inspiration cycle); and an intraoperative CBCT where the patient is in lateral decubitus position under sedation and mechanical ventilation. No surgical action was done before the CBCT acquisition, so that the patient's lung could be imaged in a fully inflated state, i.e. witouth pneumothorax.

Figure 1 shows an axial view of the images. Although the nodule is here clearly visible in both modalities, it is worth mentioning that it might not be the case, particularly for GGO nodules. While nodule visibility in intraoperative images plays an important role in image guided thoracic surgery, it is not relevant for the purposes of this work.

Segmentation

Segmentation of the lung parenchyma was achieved using an own modified version of Chest Imaging Platform, an open source library for image processing and analysis of chest CTs. [START_REF] Ross | Chest Imaging Platform: An Open-Source Library and Workstation for Quantitative Chest Imaging[END_REF] The method is based on a thresholding approach that exploits the differences in intensity of voxels inside and outside the lungs, which resulted in a segmentation containing both lungs and airways. A 3D Iterative Region Growing approach was used to segment only the primary airway branches, which were then removed from the previously obtained lung's segmentation. This allowed the separation of the two lungs, to extract the affected lung by using connected component analysis and a priori information on the position. Finally, morphological operations where applied to fill in holes and smooth irregular boundaries that were present mostly towards the mediastinum.

Given the reduced quality of the CBCT image, the aforementioned segmentation procedure was only used on the preoperative CT image. However, a similar approach using morphological eroding operations instead of airways extraction was applied to the CBCT to obtain a rough segmentation of the intraoperative lung. This latter is not accurate, and hence it is only used for the initialization of the registration workflow.

Registration

The registration process aims to account for lung deformations that occur after a change in patient position before and during a VATS intervention. A registration workflow composed of three steps is proposed: (1) initial alignment, (2) rigid registration and (3) non-rigid registration.

Initial alignment

As discussed before, preoperative and intraoperative images were taken under different configurations. Hence, the position of such images in the physical space is non-overlapping. To compensate for this misalignment, the centroids of the lung segmentations were superposed by translating the preoperative CT image. In addition, the change of orientation of the lung was assumed to be of approximately 90 degrees on the axial plane, so a rotation of this amount was also applied.

Rigid registration

Once the images were roughly aligned, an image-based rigid registration process was performed. The idea was to find a geometrical transformation that maximized image correspondence without introducing local deformations. This rigid registration step is of great importance since it affects directly the latter measurement of local deformations. Although different sources of information could be used to drive such registration (e.g. spine and ribs, airway tree, etc.), emphasis was made on structural information of the parenchyma. So, the optimization of a similarity metric estimated over the gray level information contained inside the lung mask was performed, disregarding the rest of the information on the image.

The procedure was accomplished using the multi-resolution image registration techniques implemented in the Elastix toolbox. [START_REF] Klein | elastix: A Toolbox for Intensity-Based Medical Image Registration[END_REF] Normalized Mutual Information was used as a similarity metric, with an Adaptive Stochastic Gradient Descent optimization process. For each iteration, a set of 3000 paired random points were extracted inside the lung's segmentation and were used for the computation of the similarity metric. Several image resolutions were necessary to allow the algorithm to account for large displacements.

Non-rigid registration

The last step of the registration workflow consists of maximizing image correspondence by allowing local deformations. Any non-rigid registration procedure is characterized by the similarity metric, the optimization method and the elastic transformation model. For the first two components, as for the rigid registration step, the Normalized Mutual Information similarity metric estimated on samples over the lung's parenchyma, along with an Adaptive Stochastic Gradient Descent optimizer were used. B-Splines were chosen as the elastic transformation model. Over-deformation was avoided by restricting the degrees of freedom of the transformation, i.e. reducing the amount of control points. For that, a grid spacing of 16 mm in the highest resolution was chosen, which was found empirically to be large enough to allow fine deformations but also small enough to avoid over-registration. Elastix toolbox was again used to accomplish this task, by taking the rigidly registered image as the starting point.

RESULTS

An illustration of the qualitative results obtained after applying the proposed registration workflow is presented in Figure 2. Image misalignment is represented using a complementary color approach. After the initial alignment and rigid registration steps, there is an overall good overlap of the lung's parenchyma (left). This can be particularly appreciated towards the lateral and posterior parts of the lung contour. However, important misalignments still remain towards the medial and anterior parts of the lung, since they can not be recovered without allowing local deformations. These misalignments mainly disappear after the final non-rigid registration step (middle). The lung contours are now better matched and disparities on internal structures are recovered. Before the quantification of local lung deformations, a quantitative evaluation of the registration accuracy had to be performed. For that, the Target Registration Error (TRE) between a set of 51 paired anatomical landmarks was calculated. These landmarks were manually placed by an expert thoracic surgeon both in CT and CBCT, during a single session. Instructions were given so that the spatial distribution was as homogeneous as possible, covering the whole lung. Three landmark groups could be identified: 20 landmarks on airway bifurcations, 30 landmarks on vessel bifurcations and 14 landmarks towards the periphery of the lung, where the distance to the surface is lower than 18 mm. Additionally, another landmark was placed inside the nodule. Figure 3 presents the TREs obtained before and after the non-rigid registration procedure was applied.

The spatial distribution of the anatomical landmarks can be seen in Figure 4. Emphasis was given to the set of peripheral landmarks, which was the group where the remaining TREs were the highest. It is clear from the figure that misalignments of anatomical landmarks existing after rigid registration were successfully recovered after non-rigid registration, when local deformations of the lung were taken into account.

Finally, the magnitude of the deformation field obtained after non-rigid registration is depicted in Figure 5. Two anatomical planes are used to illustrate how such deformation is distributed throughout the parenchyma. In addition, the normalized histogram of deformations present on the whole lung volume is also presented.

DISCUSSION

As for the segmentation of the lung, the proposed procedure was straightforward on the preoperative CT and produced a smooth binary mask of the lung's parenchyma. For the CBCT, however, the reduced image quality did not allow the same results. The modified segmentation procedure was then applied, and even if it could not completely recover the whole lung structure, the result was sufficient to estimate the lung's centroid required for the initialization of the registration workflow.

With respect to image registration, the TREs obtained after rigid registration only where on the order of 7 mm, with better aligned landmarks near the vessels and lung periphery (∼ 6 mm) than near the airways (∼ 8 mm). A significant error reduction was achieved after non-rigid registration. In fact, median values of TREs for all landmark groups where reduced to 0.44 mm, 0.89 mm and 1.24 mm for airways, vessels and peripheral landmarks, respectively. In particular, a TRE reduction from 7.3 mm to 0.4 mm was observed for the the landmark placed inside the nodule. It is important to note that after non-rigid registration, the largest TREs are located in the lung's periphery. Lung structures are the smallest in this regions, hence Partial Volume Effects (PVE) in the image are more important. As a result, the registration algorithm has more trouble finding correspondences, which produces larger misalignments. Overall, a median TRE of 1.01 mm was obtained for the whole set of landmarks after non-rigid registration, with a minimum of 0.35 mm and a maximum of 2.24 mm. Since the majority of the TREs are around 1.01 mm (see Figure 3), one can conclude that the resulting registration is accurate, since the errors are on the order of image spacing.

Regarding the magnitude of the deformations, experimental results showed non-negligible measurements with maximum values around 14 mm, in different areas. However, a special remark has to be made with respect to the anatomical interpretation of these deformations' locations. Although the measured deformations were here the highest near the apex and anterior parts of the lung (see Figure 5), this does not necessarily mean that these regions are the most affected by the change of configuration before the VATS intervention. In fact, the measurement of local deformations is highly sensitive to the initial alignment and rigid registration steps: a minimization of these local deformations in some regions of the lung will necessarily result in an increase in others. Hence, since there is anatomically no fixed point in the chest, deformations between the lung's areas are always relative. Other results would have been obtained with different initial geometrical transformations. However, the relative differences would remain similar.

To conclude, it is clear that important, non-uniform deformations of the lung occur between the pre-and intraoperative configurations. They are caused by a change of the patient pose, breathing mechanics and how gravity affects internal structures. Only a non-rigid registration procedure can cope with these local deformations, to obtain a perfect match between preoperative CT and intraoperative CBCT images.

CONCLUSION

This paper presents a registration workflow that was used to measure lung deformations resulting from changes between preoperative and intraoperative configurations during a VATS intervention. Registration accuracy was measured using TREs on a set of 51 paired anatomical landmarks, obtaining a median error of 7.09 mm after rigid registration, which was significantly reduced to 1.01 mm after non-rigid registration. Deformations throughout the lung where measured to be of maximum 14 mm (∼ 7 mm in average) on different regions of the lung.

Experimental results highlighted the importance of accounting for lung deformations during VATS, since it will be necessary for transforming information extracted from preoperative images to the intraoperative setting. This will be of particular importance for the localization of lung nodules, which might be not visible through intraoperative imaging when their density or size is considerably low.

Future work includes the study of lung deformations on a dataset of several cases. Also, the analysis of the deformation under different rigid initialization approaches might be of interest for the identification of deformed lung regions. Finally, special focus must be brought to the correction of deformations induced by the pneumothorax that will allow the prediction of nodule displacement during VATS, which is the long-term goal of this work.

Figure 1 .

 1 Figure 1. Slices in approximately the same transversal plane for preoperative CT (supine position) and intraoperative CBCT (lateral decubitus position). The nodule is encircled in orange. The change of configuration from preoperative to intraoperative configurations is clearly visible.

Figure 2 .

 2 Figure 2. Axial view of the lung. Left: CT-CBCT image overlap after rigid registration. Middle: CT-CBCT image overlap after non-rigid registration. Right: target intraoperative CBCT image.

Figure 3 .

 3 Figure 3. Target Registration Errors (TRE) in millimeters for different landmark groups before and after non-rigid registration.

Figure 4 .

 4 Figure 4. Spatial distribution of the whole set of landmarks. Left: Color map representation of the TREs after non-rigid registration, for all landmarks. Middle: Set of preoperative periphery landmarks (red spheres) after rigid registration compared to intraoperative periphery landmarks (blue squares). Right: Set of preoperative periphery landmarks (green spheres) after non-rigid registration compared to intraoperative periphery landmarks (blue squares).

Figure 5 .

 5 Figure 5. Left: color mapped local deformations in millimeters for coronal and axial views of the lung. Right: the normalized histogram of local deformations in millimeters for the whole lung volume.

ACKNOWLEDGMENTS

The work presented in this article was partially supported by the Région Bretagne through its Allocations de Recherche Doctorale (ARED) framework and by the French National Research Agency (ANR) through the frameworks Investissements d'Avenir Labex CAMI (ANR-11-LABX-0004) and Infrastructure d'Avenir en Biologie et Santé (ANR-11-INBS-0006).


Firas Abi-Farraj

Bernd Henze

Alexander Werner

Michael Panzirsch

Christian Ott

Máximo A Roa

Humanoid Teleoperation using Task-Relevant Haptic Feedback

Keywords: ℎ, ℎ

Robotic teleoperation is a key technology for a wide variety of fields. Teleoperating a humanoid in particular is essential as it allows the user to act remotely on an interface designed especially for humans, e.g., in a space station, or operating tools and machinery in disaster scenarios. This paper presents a 'task-relevant' haptic interface for humanoid teleoperation, which bridges the gap between the task at hand and the balance of the robot. The operator is given command over the humanoid's hands and is informed through haptic cues about the impact of her/his potential actions on the robot' stability. Moreover, a null-space autonomous controller acts in the operator's null-space to provide her/him with a wider workspace and help in the successful execution of the task. The architecture is designed to top an existing compliance controller for a torque-controlled humanoid robot. Experiments on the humanoid robot TORO are reported to demonstrate the feasibility and effectiveness of the approach.

I. INTRODUCTION

Robotic teleoperation is playing a pivotal role for allowing human presence in remote and hazardous environments such as nuclear sites, deep oceans or the outer space. Visual and haptic interfaces are usually exploited to provide the human operator with a sense of presence while giving her/him control over the robotic system. However, controlling a complex system such as a humanoid robot using human inputs alone is a challenging task that requires a lot of training and experience [START_REF] Neo | Whole-body motion generation integrating operator's intention and robot's autonomy in controlling humanoid robots[END_REF]- [START_REF] Ghalamzan | Human-in-the-loop optimisation: mixed initiative grasping for optimally facilitating post-grasp manipulative actions[END_REF]. Partial autonomy and assisted teleoperation frameworks have been proposed to reduce the mental load on the operator and improve her/his performance [START_REF] Abi-Farraj | A learning-based shared control architecture for interactive task execution[END_REF], [START_REF] Pedemonte | Visualbased shared control for remote telemanipulation with integral haptic feedback[END_REF]. Examples in this sense range from virtual fixtures [START_REF] Abbott | Haptic virtual fixtures for robot-assisted manipulation[END_REF], [START_REF] Abi-Farraj | A visualbased shared control architecture for remote telemanipulation[END_REF], which are task-dependent and require continuous human input, to more complicated shared control schemes aiming at profiting from the human's supervisory capabilities in guiding an autonomous system [START_REF] Hertkorn | Assistance for telepresence using online grasp planning[END_REF], [START_REF] Pedemonte | Visualbased shared control for remote telemanipulation with integral haptic feedback[END_REF].

On the other hand, and while many robotic systems are designed for controlled environments such as manufacturing facilities, the need for a robot capable of interacting with interfaces designed for humans is crucial. Teleoperated humanoids have been exploited as a potential solution in different applications ranging from space manipulation [START_REF] Ambrose | Robonaut: NASA's space humanoid[END_REF] to substituting humans in nuclear environments [START_REF] Iqbal | Robotics for nuclear power plants -challenges and future perspectives[END_REF], driving a lift truck [START_REF] Hasunuma | A tele-operated humanoid robot drives a lift truck[END_REF] or piloting an aircraft [START_REF] Song | Toward autonomous aircraft piloting by a humanoid robot: Hardware and control algorithm design[END_REF]. A major challenge in this regard is designing a convenient hardware-software human machine interface (HMI) capable of controlling such <firstname>.<lastname>@dlr.de highly versatile robots [START_REF] Tachi | Telexistence cockpit for humanoid robot control[END_REF], [START_REF] Yokoi | A tele-operated humanoid operator[END_REF]. The complexity of designing a teleoperation architecture for a humanoid lies in simultaneously giving the user an immersive task-oriented experience while informing her/him at the same time about the robot's balance and physical constraints. Task-relevant balance feedback is actually a topic that has not been studied thoroughly in literature, and is at the core of the work at hand. Most of the literature either provides the operator with full immersion in the task itself through haptic feedback from interactions with the environment, or gives her/him control over the posture of the humanoid while maintaining its balance. The need for bridging this gap between the robot's balance and the task at hand is indeed essential.

In [START_REF] Zucker | A general-purpose system for teleoperation of the DRC-HUBO humanoid robot[END_REF], the teleoperation system used for controlling the DRC-HUBO humanoid robot at the 2013 DARPA Robotics Challenge trials is presented. In this work, three operators were in charge of fully controlling the humanoid through a software interface for navigation and manipulation purposes. One of the operators was explicitly responsible of ensuring stability and avoiding collisions by commanding the robot through that interface.

A different approach is described in [START_REF] Neo | Whole-body motion generation integrating operator's intention and robot's autonomy in controlling humanoid robots[END_REF], where a single human operator commands a humanoid by acting on a 3 DoF (degrees of freedom) master device. The operator selects a certain point of the robot's body for manipulation, instead of simultaneously controlling all its DoFs. An autonomous controller integrates the operator's commands for producing the required whole-body motion while respecting postural stability constraints. In [START_REF] Estrelinha | Tele-operation of a humanoid robot using haptics and load sensors[END_REF], the operator is given control over the robot's posture via a PHANTOM Omni master device, over which she/he receives haptic cues, reflecting sensory data from the load cells installed at the robot's feet. The user's actions are logged and used to teach the robot how to keep balance using a learning-from-demonstration framework. A bimanual master interface using two 6 DoF master devices to give the operator command over the humanoid's legs was proposed in [START_REF] Oliveira | Bimanual haptics for humanoid robot teleoperation using ROS and V-REP[END_REF]. The operator was fed back with force cues informing her/him about the position of the Zero-Moment Point within the support polygon. The same measure was used to give the operator a sense of the robot's balance using a vibrotactile belt for providing cutaneous haptic balance feedback [START_REF] Brygo | Tele-manipulation with a humanoid robot under autonomous joint impedance regulation and vibrotactile balancing feedback[END_REF]. This cutaneous belt was replaced with a kinesthetic system applying forces on the operator's waist to inform her/him about their proximity to the edges of the support polygon in [START_REF] Ramos | A balance feedback human machine interface for humanoid teleoperation in dynamic tasks[END_REF].

The focus of the described literature is either on the manipulation task itself or on the balance/stability of the robot. In the former case, the operator is given command over the DoFs of interest while an autonomous algorithm generates the required whole-body motion and enforces stability constraints. She/he is fed with haptic information related to the task at hand, e.g., interactions with the environment, without any insight on the balance of the robot. In the latter case, the user is given the task of maintaining the robot's balance by acting on a specific HMI to control the lower body (posture, legs, ...) of the humanoid.

The work in [START_REF] Brygo | Tele-manipulation with a humanoid robot under autonomous joint impedance regulation and vibrotactile balancing feedback[END_REF] is one of few that tackles both aspects simultaneously. Here, the operator is commanding the hands of the humanoid while being informed with haptic cues of its stability through a vibrotactile belt. However, she/he is still not informed about the impact of her/his manipulative actions on stability itself. In fact, manipulation actions have a crucial impact on stabilizing/destabilizing the robot. The operator, who is given command over some of the robot's DoF, can be oblivious as to how these DoF can be employed for improving the robot's balance. To this end, we propose in this paper a balance-feedback human machine interface that closes the gap between the manipulative actions and their impact on the stability of a humanoid. A haptic interface establishes the direct association between the two by providing the operator with cues informing her/him of the impact of one on the other. These cues are meant to provide the user with potential solutions to assure a successful task completion while accounting for the different constraints of the system.

On the other hand, the null space of the operator's commands is also employed to maintain the balance of the robot through an autonomous controller acting on an underlying impedance control architecture previously presented in [START_REF] Henze | Passivity-based whole-body balancing for torque-controlled humanoid robots in multi-contact scenarios[END_REF]. In literature, a large variety of balancing controllers can be found, most of them based on inverse dynamics [START_REF] Mistry | Inverse dynamics control of floating base systems using orthogonal decomposition[END_REF]- [START_REF] Sentis | Compliant control of whole-body multi-contact behaviors in humanoid robots[END_REF]. Our approach, presented in [START_REF] Henze | Passivity-based whole-body balancing for torque-controlled humanoid robots in multi-contact scenarios[END_REF], is passivity-based and allows to stabilize the Center of Mass and the end-effectors that are not used for supporting the robot (the hands, in this context) in Cartesian space with respect to the world. This controller is then suitable as a tool for the teleoperation approach presented in this work.

The rest of the paper is divided as follows: Section II describes the modeling of the humanoid and the employed compliance controller, section III gives an account on the shared-control architecture commanding the slave's motion and the haptic interface, and the experiments and results are discussed in section IV. Finally, section V concludes the paper and discusses potential future work.

II. BACKGROUND

A. Dynamic Model

In legged humanoid robotics the use of dynamic models with a free-floating base is widespread because they feature a higher flexibility regarding contact changes compared to dynamic models with a fixed base. In general, a central body within the kinematic structure of the robot is chosen as a base link, such as the hip or the trunk. Some works also utilize the center of mass (CoM) as a base, since it represents an essential quantity for balancing. Here, we will follow the proposition of [START_REF] Henze | Passivity-based whole-body balancing for torque-controlled humanoid robots in multi-contact scenarios[END_REF] by defining a CoM frame C, which is located at the CoM and has the same orientation of the hip. Let x c ∈ R 3 and R c ∈ SO(3) denote the position and orientation of the frame C with respect to the world frame W. The corresponding translational and rotational velocities are ẋc and ω c , respectively. Based on the n joint angles q ∈ R n and v c = ( ẋT c , ω T c ) T , the dynamics of the humanoid robot is given by

M vc q ν +C v c q ν + mg 0 0 g = 0 τ + τ ext . (1)
Herein, M ∈ R (6+n)×(6+n) and C ∈ R (6+n)×(6+n) denote the inertia and Coriolis/centrifugal matrix, respectively. The gravitational torques are given by g ∈ R 6+n with m denoting the overall mass of the robot and g 0 ∈ R 6 the gravitational acceleration 1 . The joint torques are given by τ ∈ R n . The influence of external wrenches acting on the robot is taken into account by the generalized torque vector τ ext ∈ R 6+n . Let us divide the Ψ end-effectors into two subgroups [START_REF] Henze | Passivity-based whole-body balancing for torque-controlled humanoid robots in multi-contact scenarios[END_REF]: The first one is referred to as "balancing end-effector" (bal) and contains the ψ end-effectors that are used by the robot to support itself (usually the feet). The remaining end effectors are called "interaction end-effectors" (int), as they are still free to be used in a manipulation or interaction task (usually the hands). Based on this definition, the Cartesian velocities of the end effectors v ∈ R 6Ψ are given by

v = v bal v int = J bal J int ν = J ν. (2)
with the Jacobian matrix

J ∈ R 6Ψ×(6+n) and v bal ∈ R 6ψ , v int ∈ R 6(Ψ-ψ) , J bal ∈ R 6ψ×(6+n) , J int ∈ R 6(Ψ-ψ)×(6+n) .
In the case where all external disturbances act solely at the end-effectors, τ ext simplifies to

τ ext = J T F ext (3)

B. Underlying Compliance Controller

This section gives a brief recapitulation of our balancing controller presented in [START_REF] Henze | Passivity-based whole-body balancing for torque-controlled humanoid robots in multi-contact scenarios[END_REF]. The controller stabilizes the CoM by a Cartesian compliance, which applies a wrench F c ∈ R 6 at the CoM frame C . Each one of the interaction end-effectors is stabilized by another Cartesian compliance, with the resulting wrenches stacked into F int ∈ R 6(Ψ-ψ) . In order to support the robot, the control algorithm computes a suitable set of balancing wrenches F bal ∈ R 6Ψ by solving the following quadratic optimization problem

F opt bal = argmin F bal F bal -F def bal T Q F bal -F def bal ( 4
)
with respect to

Ad T bal F bal + Ad T int F int = mg 0 -F c (5)
and

f k,z ≥ f min k,z , δ k ∈ S k , f k,x/y ≤ μk f k,z . (6)
The cost function (4) minimizes the deviation of F bal from a default wrench distribution F def bal considering the positive definite weighting matrix Q ∈ R 6Ψ×6Ψ . The default distribution F def bal is a tuning parameter, which can be set to half of the weight of the robot if the stance is symmetric. The equality constraint (5) represents the underactuation of the base by demanding that the influence of all commanded end-effector wrenches (F bal , F int ) on the CoM must sum up to the compliance wrench F c plus gravity. For this, the Jacobian matrices

J bal = [ Adbal Jbal ] and J int = [ Adint Jint ] are partitioned into Ad bal ∈ R 6ψ×6 , Ad int ∈ R 6(Ψ-ψ)×6 , J bal ∈ R 6ψ×n
, and J int ∈ R 6(Ψ-ψ)×n . The first two are the stacked adjoint matrices of each end-effector, and relate a motion of the CoM frame C with a motion of the endeffectors. The matrices J bal and J int describe the influence of a joint motion on the end-effectors. The inequality constraints ( 6) represent the contact model to which F bal is subjected to in order to account for unilaterality, the position of the Center of Pressure (CoP), and for friction of the balancing contacts. For each wrench within F bal , the force perpendicular to the contact surface S is bounded from below by the minimum contact force f min k,z in order to prevent the end-effector from lifting off 2 . Slippage is prevented by constraining the tangential forces f k,x/y to the friction cone given by μk . The CoP δ k is restricted to lie inside the contact surface S to prevent the end-effector from tilting. From the contact model ( 6) and a given contact configuration, one can compute an equivalent support polygon, as shown in [START_REF] Caron | Multi-contact walking pattern generation based on model preview control of 3D COM accelerations[END_REF]. In order to achieve a static and stable balancing, the CoM projected to the ground floor must stay within the support polygon as well. This allows us to use the support polygon as a stability criterion in Sec. III-C. After computing a suitable wrench distribution F opt bal , the end effector wrenches are mapped to joint space via

τ = -J T bal J T int F opt bal F int . (7)

III. SHARED-CONTROL ARCHITECTURE

A. Whole-Body Motion

While it is habitual to have the high-level 'manipulation' controller incorporated into the compliance controller itself, a different approach is opted for in this work. Indeed, the compliance controller described in the previous section is treated as a 'packaged' controller and topped with the needed higher level position/velocity whole-body motion controller, as shown in Fig. 2. This approach is remunerative as it allows the user to design a controller for a force/torque compliancecontrolled robot from a pre-designed and tuned compliance controller without having to deal with its specificities.

Going back to Fig. 2, the shared-control teleoperation architecture is the upper gray zone. The user acts on an HMI sending velocity commands (v m,r , v m,l ) through the right and left arms of the master device while receiving haptic cues (F m,r , F m,l ) from potentials informing about the proximity to the edge of the support polygon and other physical constraints of the system. As the user gets closer to hitting these constraints, her/his commands are gradually scaled down to zero in order to ensure the stability of the system. On the other hand, the mentioned potentials are also used to generate a CoM motion v c in the null-space of the motion of the hands, which are commanded by the user, and the feet (stable on the ground) in order to ensure that the system remains as far as possible from the constraints. The resulting velocities are integrated to generate new positions, which are then fed to the underlying compliance controller.

B. Master Side

As described before, and following the classical bilateral force-feedback teleoperation framework, we assume the presence of a master device consisting of a right and a left arm, through which the user sends velocity commands and receives force-feedback. The two master arms are modeled as generic (gravity pre-compensated) mechanical systems

M m (x m )ẍ m + C m (x m , ẋm ) ẋm = F m,r F m,l + F h (8)
where x m ∈ R 12 is the device cartesian configuration vector containing the pose of both the right and left master arms, M m (x m ) ∈ R 12×12 is the positive-definite symmetric mass matrix, C m (x m , ẋm ) ẋm ∈ R 12×12 accounts for Coriolis/centrifugal terms, F h ∈ R 12 account for the forces applied by the human operator, and F m,r , F m,l ∈ R 6 are the control forces on the right and the left arm, respectively. On the other hand, the individual linear and angular velocities of each of the master's right and left arms' end-effectors are denoted by v m,r ∈ R 6 and v m,l ∈ R 6 respectively. The operator is given control over the hands of the robot through a direct Cartesian coupling. As the workspace of a human operator and the master arms is close to that of the humanoid, no scaling is implemented. However, the user can always "clutch" to move the master arms to a more convenient position without moving the slave. A velocityvelocity master-slave coupling is employed such that

v h,r = (1 -α) Rk v m,r v h,l = (1 -α) Rk v m,l (9)
where v h,r and v h,l ∈ R 6 are the desired velocities of the right and left hands of the slave respectively, and Rk is the rotation map between the master and the slave velocities. α : 0 → 1 is an activation function that goes from 0 to 1 as the system approaches any of the different constraints, with 0 being away from the constraint and 1 being exactly at it.

C. Balancing and Physical Constraints

In addition to executing the operator's commands, it is crucial to maintain the stability of the humanoid by ensuring its CoM remains within the support polygon 3 , as shown in Fig. 3. To this end, we define a cost function H b : R n → R as a measure of the 'balance' of the humanoid, which goes to infinity as the user approaches the edges of the support polygon such that

H b (p c ) = i h i,b (p c ), (10)
where p c ∈ R 3 is the position of the center of mass, and h i,b is the cost function attached to each edge i of the support polygon and defined by 3 Note that this condition is only valid for evaluating static stability. In the dynamic case, a condition could be used based e.g. on the capture point [START_REF] Pratt | Capture point: A step toward humanoid push recovery[END_REF].

h i,b (p c ) =    µ tan π 2 1 -di d l,i 2 d i < d l,i , 0 otherwise (11)
with µ being a regulation gain (µ = 1 in this work), d i is the horizontal distance from the center of mass to each edge i of the support polygon (which is known from the fixed posture of the feet), and d l,i is the distance from the center of the support polygon to its i-th edge.

x The time derivative of H b can then be defined by

Ḣb = ∂H b ∂p c ṗc = ∂H b ∂p c J p c ν. (12)
where J p c = [I 0] ∈ R 3×6+n is the jacobian mapping the state velocities ν to ṗc . While the physical constraints of a humanoid robot range from workspace limits to joint limits, singularities and torque constraints, a representative workspace constraint was considered in this work consisting of the distance between the hands of the robot and its center of mass. Constraining this distance ensures that the robot maintains a "safe" posture and avoids a "full stretch" configuration that could push the robot to singularities and joint limits. While this constraint does not cover the various limitations of a humanoid, it is meant as a proof of concept and a representative constraint. The exploration of further constraints is a future work.

The constraint hereby considered is then to maintain the distance between each of the hands of the humanoid and its CoM within a predefined threshold d th . The corresponding cost function H d accounting for this physical constraint is defined as [START_REF] Yokoi | A tele-operated humanoid operator[END_REF] where the cost function of each hand h c,x p c , p h,x is

H d p c , p h,r , p h,l = h c,r p c , p h,r + h c,l p c , p h,l
h c,x =    ρ tan π 2 1 - d h,x -ds d th -ds 2 d s < d i < d l , 0 otherwise .
(14) d h,x in the above equation is the distance from each hand x to the center of mass, while d s is the minimum distance after which the potential starts, and ρ is a positive gain (ρ = 1 in this work). The derivative of the potential H d is then

Ḣd = ∂H d ∂p c ∂H d ∂p h,r ∂H d ∂p h,l   ṗc ṗh,r ṗh,l   = ∂H d ∂p c ∂H d ∂p h,r ∂H d ∂p h,l   J p c J p h,r J p h,l   ν, (15)
where J p h,r and J p h,l are the jacobians linking the linear velocities of the right and left hand to the state velocity ν, respectively.

The desired state velocity ν des to minimize the potentials and make sure the system constraints are respected can then be defined as the negative of the transpose of the gradient of potentials H b and H d (eq. ( 12) and eq. ( 15)) such that

ν des = -   ∂H d ∂p c ∂H d ∂p h,r ∂H d ∂p h,l   J p c J p h,r J p h,l     T - ∂H b ∂p c J p c T (16)

D. CoM Control Law

The motion of the center of mass is designed to follow the commands of the human operator while fixing the feet on the ground as a primary task, and maintaining the balance of the robot in their null-space. We commence now by explaining the autonomous balancing algorithm governing this behavior.

The robot's state velocity vector is mapped to the velocities of the hands and the feet (which are stationary in the case at hand) through a jacobian

J     v h,r v h,l v f,r v f,l     = J v c q . ( 17
)
The impact of the velocity inputs from the human operator on the velocity of the CoM, dubbed as its primary task, can then be retrieved by inverting the previous equation such that

v c = Ad †     v h,r v h,l v f,r v f,l     ( 18
)
with Ad † being the Moore pseudo-inverse of

Ad = [ Ad T bal Ad T int ]
T , which is the upper part of Jacobian J corresponding to the velocity of the center of mass v c . The null-space balancing motion is then added to [START_REF] Brygo | Tele-manipulation with a humanoid robot under autonomous joint impedance regulation and vibrotactile balancing feedback[END_REF], and the resulting full control law becomes

v c = Ad †     v h,r v h,l v f,r v f,l     + N ν des ( 19
)
where N = (I -J J † ) is the null-space projector.

E. Haptic Feedback

In case of proximity to a constraint, force-feedback is given to the operator to guide her/him away from the undesired configuration. Moreover, if the operator keeps pushing towards the constraint, her/his velocity commands are tuned down to prevent the system from hitting the constraint.

A proximity function α = max(α i (d i )) is employed for that purpose, where α i : 0 → 1 is the proximity function associated to each constraint i. The constraints in this case include the different edges of the support polygon and the distances from the hands to the CoM. α i (d i ) increases linearly from 0 to 1 as a function of the distance to the corresponding constraint d i , such that it is 0 before a predefined threshold and reaches 1 at the constraint

α i (d i ) : 0 → 1 d i : d p,i → d l,i . ( 20
)
with d p,i being a maximum threshold after which α i becomes active, and d l,i is the limit of the constraint.

On the other hand, the desired motion direction of the hands [ṽ h,r , ṽh,l ] that ensures the minimization of the potentials H b and H d can be retrieved as

ṽh,r ṽh,l = - J h,r J h,l N bal ν des , (21)
where

N bal = (I -J bal J † bal )
is the null-space projector ensuring that the generated motion does not impact the position of the feet. J bal is the jacobian mapping the state velocity ν to the velocity of the feet

[v f,r , v f,l ].
The input velocities of the human operator are then tuned down if the robot is in the proximity of a constraint, i.e. α > 0, and the direction commanded by the operator is opposing the desired motion directions [ṽ h,r , ṽh,l ]

v h,r = (1 -α) Rk v m,r if ṽT h,r Rk v m,r < 0 v h,l = (1 -α) Rk v m,l if ṽT h,l Rk v m,l < 0 . (22)
The forces received by the operator on the right and left arms are defined by

F m,r = F max α ṽh,r /||ṽ h,r || + B v m,r F m,l = F max α ṽh,l /||ṽ h,l || + B v m,l (23)
where F max is a design parameter defining the maximum force expected from the operator, and B ∈ R 6×6 is a damping matrix. This design of the haptic interface is unique in that it is not a pure resistive force stopping the user as she/he approaches a constraint or another system limit. On the contrary, the described force cues are an active guidance providing the user with several solutions for avoiding the system constraints over the 12 DoF she/he is commanding.

Note that the behavior of the force feedback in eq. ( 23) is akin to a spring-damper system, as α is directly proportional to the distance between the CoM and the support polygon boundaries. The potentials H b and H d do go to infinity in the proximity of the edges of the support polygon, but the haptic feedback given to the operator adopts the direction of the gradient of these potentials (but not its magnitude). While this is not needed for passivity, as the potential can be proven passive in itself, it allows for a better design of the haptic interface, since any haptic device has limited capabilities for force generation and the design of the haptic interface must account for these limitations.

IV. EXPERIMENTS AND RESULTS

This section describes the experiments conducted to validate the described approach, and discusses the obtained results.

Our experimental setup uses HUG on the master side. HUG is a haptic user interface system composed by two light-weight robotic arms with a workspace of 1.1 m and a nominal payload of 20 kg each [START_REF] Sagardia | A platform for bimanual virtual assembly training with haptic feedback in large multi-object environments[END_REF]. On the slave side, we use TORO, a humanoid robot developed at DLR [START_REF] Englsberger | Overview of the torque-controlled humanoid robot TORO[END_REF]. It has a height of 1.74 m, a weight of 76.4 kg, and 25 joints in total (not counting the neck and the hands). In the presented experiments, the feet of TORO are in contact with the ground floor in order to support the robot (balancing end-effectors). The hands (interaction end-effectors) were directly commanded by the operator via HUG. The joints in the neck and hands were not used. Due to the feet contacts (2×6 = 12 DoF) and the user input for the hands (2×6 = 12 DoF), the robot can only move 7 DoF in order to maintain balance (6 DoF attached to the floating-base (frame C) plus 1 remaining DoF).

In order to test the behavior of the different parts of the system, we distinguish two major components: the nullspace autonomous balancer acting on the CoM to keep the system as far as possible from constraints, and the haptic guidance informing the user of the proximity to any of these constraints and the possible directions to avoid them. Three experiments were performed while activating or deactivating these two components, to asses their impact on the general behaviour of the system. The video attached to this paper provides further insight on the performed experiments.

A. Experiment I

In the first experiment, the haptic guidance and the nullspace balancer were both deactivated, and the user was given unrestricted control over the hands of the humanoid.

The user was asked to reach as far as possible with the hands of the humanoid along the x-axis (which is pointing forward). Fig. 4 reports on the obtained results. Fig. 4a and Fig. 4b show the x-position of the hands and the CoM respectively throughout the experiment, while Fig. 4c shows the variation of the potentials H b and H d . Note that this is just a visualization of the potentials, as they were not active in this experiment. On the other hand, in Fig. 4b the horizontal dotted line represents the edge of the support polygon along the positive x-direction. The center of the support polygon was at 1 cm from the world frame in this scenario, with its edges at -3 cm and 7 cm along the x-axis. Observing the figures, the maximum reach of the hands of the humanoid before hitting the support polygon's edge was 52.7 cm for the right arm (the vertical dotted line denotes the moment at which the CoM crossed the edge of the support polygon). In fact, the humanoid then started to tip over, losing contact on its right foot, and the experiment was called to a halt. Note that the potential H b increased significantly as the edge of the support polygon was closer.

The results of this experiment are also depicted in Fig. 7, where Fig. 7a shows the trajectories followed by the right hand (in blue) and the CoM (in red), and Fig. 7d depicts the final posture of the robot after the right foot lost grip.

B. Experiment II

Following on the previous experiment, haptic guidance was activated along with the restriction applied on the operator's commands when approaching a constraint (by exploiting the proximity measure α, check eq. ( 22)). The nullspace autonomous balancer was, however, still not active. The same experiment was repeated and the user was asked again to reach as far as possible. The results are reported in Fig. 5. Fig. 5b plots the forces fed to the operator (along the x-axis) on the end-effectors of the master device. On the other hand, Fig. 5a and Fig. 5c show the x-position of the slave hands (right and left) and the center of mass, respectively, as the experiment proceeds. Finally Fig. 5d depicts the two potential functions describing the stability of the robot and its physical constraints.

An interesting behavior is observed in this experiment. As the user moved the hands of the robot forward, the CoM started approaching the edge of the support polygon, and this was reflected as an increase in the cost functions and an active force on both hands along the negative x direction (see t=10 sec). The user then reactively started moving the left hand backwards, following the haptic guidance. This decreased the cost function, allowing him to further push the right hand forward. In fact, this strategy, inspired by the informative haptic guidance, helped the user to reach 34.9 cm farther than the first experiment (t=45 sec), without impacting the stability of the robot. The followed trajectories and final posture of the robot for experiment II are reported in Fig. 7b and Fig. 7e.

C. Experiment III

The null-space autonomous balancer was activated for the third experiment, in addition to having the haptic guidance active as well. The experiment followed the same procedure as before, and the results are reported in Fig. 6.

An interesting figure to observe in this experiment is Fig. 6c, which shows the evolution of the x-position of the CoM. In contrast to the previous two experiments, the CoM remained close to the center of the support polygon for a significant portion of the experiment. In fact, both the right and left hands were close to a 60 cm reach before the CoM started moving forward. It was then pushed forward as the hands were being moved farther ahead, thus driving the robot to be more stretched and increasing the potential H d , as can be observed in Fig. 6d. The resulting maximum reach was 91.4 cm, a bit higher than in experiment II. However, the user was able to push both hands forward instead of one, thus allowing for more manipulation options. On the other hand, Fig. 6b shows that the forces received throughout this experiment were negligible, except for t= [33,48] s when the hands were stretched at more than 85 cm. This is important in that autonomy, while accounting for all the operator's commands, was able to successfully command the null-space and provide the operator with the maximum workspace possible with minimal disturbance.

The followed trajectories and final posture of the robot in this experiment are reported in Fig. 7c, and Fig. 7f.

V. CONCLUSIONS

This paper presented a shared-control architecture for teleoperating a torque-controlled humanoid robot. An operator was given control over the hands of the humanoid, while the whole-body motion was governed by a null-space balancer acting in the null-space of the operator's commands. A novel approach for providing haptic feedback is introduced, where the user is fed with high-level informative haptic cues informing her/him about the impact of her/his potential actions on the robot's balance. This approach bridges the gap between the task itself and the different constraints of the system, thus allowing the operator to adapt her/his approach for a successful task execution within the constraints of the system. On the other hand, the authors are looking forward to explore different constraints of the system both in Cartesian space (collision avoidance) and in joint space (joint limits, singularities, ...), and study the possibility of incorporating them in the architecture in the future. Another potential future work is the consideration of variable support polygons, which change depending on the (changing) posture/configuration of the robot. Finally, incorporating dynamics into the approach to provide the operator with a meaningful instantaneous guidance during a dynamic behavior, like walking, is another future challenge.

Fig. 1 .

 1 Fig. 1. The Experimental Setup. The haptic user interface system HUG (right) is used as master device, while the torque-controlled humanoid robot TORO (left) works as slave device.

Fig. 2 .

 2 Fig.2. A scheme of the proposed shared-control framework. The lower block represents the underlying compliance controller, while the upper block is the high-level teleoperation framework.

Fig. 3 .

 3 Fig. 3. An illustration of the support polygon for a robot standing on horizontal ground.

Fig. 4 .

 4 Fig. 4. Experiment I: No haptic guidance is provided, and the null space balancer is deactivated. (a) and (b) show the position of the hands and CoM along the x-axis, respectively, while (c) shows the system potentials H b and H d .

Fig. 5 .

 5 Fig. 5. Experiment II: Haptic guidance is provided to the operator, but the null-space balancer is deactivated. (b) shows the force cues fed to the operator along the x-axis on the right and left hand, respectively, while (a) and (c) report on the x-position of the hands and the CoM. (d) plots the potentials H b and H d .

Fig. 6 .

 6 Fig. 6. Experiment III: Haptic guidance is provided to the operator, and the null-space balancer is active. (b) shows the force cues fed to the operator along the x-axis on the right and left hand, respectively, while (a) and (c) show the x-position of the hands and the CoM. (d) depicts the potentials H b and H d .

Fig. 7 .

 7 Fig. 7. Comparison of the conducted experiments: The top row shows the trajectory of the right hand and the CoM in the saggital plane of the robot. The bottom row shows the final pose of the robot with the largest achievable reach of the hands. (a) and (d) report on experiment I, (b) and (e) on experiment II, and (c) and (f) on experiment III

The author is with CNRS at Irisa and Inria Rennes Bretagne Atlantique, Campus de Beaulieu, 35042 Rennes Cedex, France.firas.abi-farraj@irisa.fr

[START_REF] Hertkorn | Assistance for telepresence using online grasp planning[END_REF] The authors are with the German Aerospace Center (DLR), Institute of Robotics and Mechatronics, 82234 Wessling, Germany.

Note that g 0 is six-dimensional by containing also the rotational DoFs.

For the conducted experiments (see Sec. IV), f min k,z was set to 50 N in order to account for joint friction.

ACKNOWLEDGMENT

This work has received funding from the European Union Horizon 2020 Research and Innovation programme under grant agreement No. 645582 (RoMaNS), and H2020-ICT-645097 (COMANOID).