Gautier Laurent
email: gautier.laurent@univ-orleans.fr

Guillaume Caumon

Mark Jessell

Interactive Editing of 3D Geological Structures and Tectonic History Sketching via a Rigid Element Method

Keywords:

Highlights

A Displacement interpolation process 18

B Linear system computation details 18

1 Introduction 3D models of geological structures are a fundamental ingredient of natural resource characterisation. Indeed, accurate modelling of geological structures is crucial for the proper characterisation of the geometry and the connectivity of rock volumes. Such structural models are generally obtained by a static and deterministic procedure, which relies on geometrical compatibility with geophysical data and well interpretations [START_REF] Mallet | Geomodeling[END_REF][START_REF] Caumon | Surface-based 3D modeling of geological structures[END_REF]. Kinematical and mechanical consistency of the model can be achieved when dense, high quality data are present and through appropriate decisions during modelling, but is not enforced by the modelling algorithms themselves.

Restoration of the geological model is a way to check these aspects a posteriori, and also to provide critical insights into the tectonic history, but computation time and modelling effort may be prohibitive [START_REF] Moretti | Working in complex areas: New restoration workflow based on quality control, 2D and 3D restorations[END_REF][START_REF] Durand-Riard | Balanced restoration of geological volumes with relaxed meshing constraints[END_REF], 2013]. When inconsistencies are detected, the model may have to be edited or rebuilt, which may lead to significant additional modelling efforts. Taking the deformation history into account at an early stage of the 3D structural modelling appears to be a natural way to produce geologically consistent models, which should reduce the need for quality control.

When building a model, several tectonic scenarios may have to be rapidly tested, which calls for tools that are simple, efficient and robust enough to be used in an automated process, while allowing for rapid user feedback and editing. In this paper, we introduce new specifications for a general dynamic deformation tool to let geologists sketch and rapidly test deformation scenarios on 3D geomodels, for example to build synthetic 3D geometries based on a succes-sion of events (as in [START_REF] Jessell | Structural geophysics: integrated structural and geophysical modelling[END_REF], restore a 3D geomodel or edit a 3D structural interpretation (Section 2). These specifications are based on two key ideas:

1. The scale at which the deformation has to be modelled is generally larger than the scale of the model. Depending on the application, the scale and resolution of a geological model may depend on lithological and stratigraphic heterogeneities, and on local geometrical peculiarities (e.g. thin stratigraphic layers, acute angular contacts, local features such as persistence of paleo-relief). The deformations that are taken into account in the modelling generally impact the global shape of the geological model while preserving these detailed features. Therefore the wavelength of the deformation may be considerably larger than the size of the mesh, which is used to describe the geometry of the model. 2. Because the full deformation path is generally poorly constrained, the interactivity of the deformation tool may be more important than its physical accuracy.

Certainly, accurate simulation of rock deformations are required in certain applications (e.g. in studies aiming at improving theoretical understanding of tectonic principles), however when few data are available while building quantitative structural models, the benefit of an accurate simulation would be limited by the associated uncertainties. On the contrary, in an exploration context, rapidly testing several deformation scenarios with approximate simulations could be a useful step before performing more precise, but more expensive, simulations. To meet these specifications, we propose to use a dynamical deformation model (Section 3), first presented in the Computer Graphics community [START_REF] Botsch | Adaptive space deformations based on rigid cells[END_REF]. It consists of a spatial discretisation with rigid polyhedrons that are linked by a non-linear cost function similar to the elastic mechanical energy. This cost function measures the difference in displacement between neighbouring elements. Boundary conditions are set on certain elements in order to honour typical subsurface data. The final deformation is obtained by optimising the displacements with respect to the cost function. The displacement field is finally applied to the geological objects embedded into the rigid elements.

In Sections 4 and 5, we show a few examples of the capabilities of this deformation algorithm to build, edit and restore geological structures. We also illustrate how this approach can be applied to a complete forward modelling exercise. With this approach, 3D structural models would then be produced by successively applying the effects of kinematic events on an assumed initial geometry until the observations of the current structures are honoured [START_REF] Jessell | Structural geophysics: integrated structural and geophysical modelling[END_REF]. The underlying tectonic history of resulting models is thus explicitly controlled by the interpreter and can be used to study structural uncertainties.

Deformation algorithms in geomodelling

In this section, we analyse the different uses of deformation algorithms for modelling geological structures and point out the need for a new range of algorithms which reinforce robustness and interactivity compared to physical accuracy. We derive new specifications from these observations and briefly review existing algorithms.

The uses of deformation computation for structural modelling

Deformation tools are, or may be, used at different stages throughout the construction of a 3D structural model: Restoration: this process aims at assessing the compatibility of a given structural model by unfolding and unfaulting the geological structures, which also provides precious insights into their possible paleo-geometry. Historically, the first approaches to restoration aimed at restoring cross-sections based on geometrical criteria [START_REF] Chamberlin | The Appalachian folds of Central Pennsylvania[END_REF][START_REF] Dahlstrom | Balanced cross sections[END_REF]. Since then, this method has been extensively developed and adapted to 3D geological structures through a variety of approaches: continuum theory [START_REF] Mallet | Geomodeling[END_REF][START_REF] Massot | Implémentation de méthodes de restauration équilibrée 3D[END_REF], surface unfolding [START_REF] Gratier | Restoration and balance of a folded and faulted surface by best-fitting of finite elements: principle and applications[END_REF][START_REF] Williams | Restoration and balance of complex folded and faulted rock volumes: flexural flattening, jigsaw fitting and decompaction in three dimensions[END_REF][START_REF] Egan | Three-dimensional modelling and visualisation in structural geology: New techniques for the restoration and balancing of volumes[END_REF][START_REF] Rouby | 3-D restoration of complexly folded and faulted surfaces using multiple unfolding mechanism[END_REF][START_REF] Griffiths | A new technique for 3-D flexural-slip restoration[END_REF], mechanical methods [START_REF] Maerten | Chronologic modeling of faulted and fractured reservoirs using geomechanically based restoration: Technique and industry applications[END_REF][START_REF] Moretti | Kine3D: A new 3D restoration method based on a mixed approach linking geometry and geomechanics[END_REF][START_REF] Muron | Méthodes numériques 3D de restauration des structures géologiques faillées[END_REF]. The application of 3D kinematic methods is relatively easy but relies on assumptions about the deformation mode (e.g. simple shear or flexural slip) and relies on global projections to approximate volumetric deformations. In contrast, geomechanical restoration is much more flexible but its application is more difficult due to preliminary conformal meshing and the possibly large number of mechanical parameters to be defined. In this paper, we will show how a light volumetric deformation tool can be used to fill the gap between these methods (Section 5.3).

Forward History Modelling: deformation algorithms may also be used to go forward in time. Starting from a hypothetical initial geometry, several kinematic deformation events can be progressively applied to pre-existing geometries. This method allows us to build very complex 3D structures with historical and geological control [START_REF] Jessell | Structural geophysics: integrated structural and geophysical modelling[END_REF][START_REF] Caumon | Current bottlenecks in geomodeling workflows and ways forward[END_REF]. Advanced methods for physically-based forward modelling of large deformations also exist, for instance Arbitrary Lagrangian-Eulerian schemes [START_REF] Allken | Three-dimensional numerical modeling of upper crustal extensional systems[END_REF][START_REF] Loiselet | Subducting slabs: Jellyfishes in the Earth's mantle[END_REF].

Local Editing: deformation may be applied at a more local scale to modify the geometry of a model, for example for integrating new data or interpretations. It generally involves locally rebuilding a part of the model after modifying some data or constraints [START_REF] Grosse | Mise en cohérence automatique d'un modèle géologique 3D[END_REF][START_REF] Caumon | Building and editing a Sealed Geological Model[END_REF][START_REF] Tertois | Editing faults within tetrahedral volume models in real time[END_REF].

Uncertainty Space Exploration: the deformation applied to a model may be used to explore the related uncertainties, for example by adding faults and consistently deforming surrounding objects [START_REF] Georgsen | Fault displacement modelling using 3D vector fields[END_REF][START_REF] Cherpeau | Method for stochastic inverse modeling of fault geometry and connectivity using flow data[END_REF][START_REF] Laurent | A parametric method to model 3D displacements around faults with volumetric vector fields[END_REF]. A history matching process, ensuring the consistency with production data, may also be applied to reduce the geological uncertainties, which implies editing the model in order to minimise the production misfit. This process generally focuses on the rock properties, but it can also include structural aspects [START_REF] Cherpeau | Method for stochastic inverse modeling of fault geometry and connectivity using flow data[END_REF][START_REF] Seiler | Structural surface uncertainty modeling and updating using the ensemble Kalman filter[END_REF]. These validation steps come at a high computational cost and require several iterations of the modelling procedure. As a result of the uncertainties associated with the available data, their interpretation and the geological processes, several possible models should be produced to properly explore the domain of uncertainties and the relevant deformation scenarios. De-formation algorithms can then be used to ease the construction of new possible models by modifying existing ones.

Specifications for a light geological deformation algorithm

Deformation algorithms are principally used in geology to run precise geomechanical simulations, either for restoration [START_REF] Durand-Riard | Balanced restoration of geological volumes with relaxed meshing constraints[END_REF][START_REF] Maerten | Chronologic modeling of faulted and fractured reservoirs using geomechanically based restoration: Technique and industry applications[END_REF][START_REF] Moretti | Kine3D: A new 3D restoration method based on a mixed approach linking geometry and geomechanics[END_REF] or forward modelling [START_REF] Burbidge | Numerical models of the evolution of accretionary wedges and fold-and-thrust belts using the distinct-element method[END_REF][START_REF] Schmid | Evolution of large amplitude 3D fold patterns: A FEM study[END_REF][START_REF] Loiselet | Subducting slabs: Jellyfishes in the Earth's mantle[END_REF][START_REF] Allken | Three-dimensional numerical modeling of upper crustal extensional systems[END_REF]. These simulations rely on information about the mechanical material and detailed geometry, which may be hypothetical or unavailable at the early stage of structural model construction. In fact, for most of the applications listed in the previous section, the physical accuracy may not be the most important aspect due to the lack of observations. This observation leads us to define the following specifications for a light generic deformation tool for geological editing and deformation scenario testing, with a relaxed constraint on physical accuracy. The following criteria are presented in order of importance:

Interactivity: the algorithm has to be fast enough to allow interactive deformation. For example, it should provide a result in at most a couple of seconds when directly manipulated by the user for editing a geological model. This way the user could interactively modify the boundary conditions while the deformation is recomputed. For more complex deformation, it has to allow us to explore different deformation scenarios in an acceptable time by iterative algorithms such as response surfaces and stochastic approaches. This would require us to compute hundreds of deformation scenario in a few minutes.

Adapted Scale: the deformation has to be computed at a scale relevant of its dimension. The deformation algorithm should not be affected by small-scale heterogeneities and geometrical details unless they lead to localisation phenomena.

Robustness: the deformation computation has to stay interactive and avoid breaking down even under extreme deformations. In particular, constitutive elements used for computation (typically triangles or tetrahedrons) should never degenerate during the process, i.e. their geometry and neighbourhood relationships should be maintained allowing the optimisation process to continue.

Topological Self-Consistency: the deformation has to affect all the structural objects contained in the model without perturbing their structural relationships. For discontinuities like faults and unconformities, the requirement to preserve these characteristics depends on the application context.

Parsimony: the 3D model has to be manipulated through a small set of parameters, both to simplify user interaction and to limit the problem of empty space for inverse approaches.

Physical Consistency: the deformations produced have to be kinematically consistent, i.e. honour possible natural displacement fields for both small and large deformations, while minimising deformations, avoiding self-intersections and local collapses or volumetric inversions. Ideally, the implied mechanical behaviour should also be consistent with natural geological materials, a condition that may be relaxed because of the large uncertainties about initial geometries and materials involved in the deformation. A first order compatibility may be enough, in particular because relaxing this constraint allows us to potentially gain in interactivity and robustness.

Extensibility: the algorithm's formalism has to remain adaptable to geological discontinuities, such as faults, and the specific mechanical behaviour of rocks (e.g. transverse isotropic).

Reversibility: the algorithm should be able to reverse the deformations it produces. This does not imply that the deformation has to be physically reversible, but the algorithm should be able to go backward as long as all the parameters of the deformation are known. This is necessary to interactively add and remove deformation components.

Overview of physically-based deformation algorithms

A large variety of deformation algorithms exists in the literature, mainly because of the diversity of applications and the variability of deformation behaviour. Here, we focus on mechanical models that are based on some physical considerations allowing us to control the kinematical and mechanical aspects of the deformation.

In the context of geomodelling, the most common deformation methods are the Finite Element Method (FEM), mainly for restoration purposes [START_REF] Durand-Riard | Balanced restoration of geological volumes with relaxed meshing constraints[END_REF][START_REF] Maerten | Chronologic modeling of faulted and fractured reservoirs using geomechanically based restoration: Technique and industry applications[END_REF][START_REF] Moretti | Kine3D: A new 3D restoration method based on a mixed approach linking geometry and geomechanics[END_REF][START_REF] Muron | Méthodes numériques 3D de restauration des structures géologiques faillées[END_REF][START_REF] Schmid | Evolution of large amplitude 3D fold patterns: A FEM study[END_REF][START_REF] Siavelis | Modélisation numérique X-FEM de grands glissements avec frottement le long d'un réseau de discontinuités[END_REF][START_REF] Vidal-Royo | Multiple mechanisms driving detachment folding as deduced from 3D reconstruction and geomechanical restoration: the Pico del Águila anticline (External Sierras, Southern Pyrénées)[END_REF], the Discrete Element Method (DEM), mainly for studying rock deformation and tectonics [START_REF] Burbidge | Numerical models of the evolution of accretionary wedges and fold-and-thrust belts using the distinct-element method[END_REF][START_REF] Scholtès | A DEM model for soft and hard rocks: Role of grain interlocking on strength[END_REF][START_REF] Vidal-Royo | The roles of complex mechanical stratigraphy and syn-kinematic sedimentation in fold development: insights from discrete-element modelling and application to the Pico del Águila anticline (External Sierras, Southern Pyrenees)[END_REF], and Arbitrary Lagrangian-Eulerian Methods (ALE). They generally aim at producing an accurate simulation of physical rock behaviour.

In a different context, the Computer Graphics community proposes a very different view of deformation computation: it has to be accurate enough to be realistic but fast enough to be interactive. In other words, the requirement of interactivity and robustness may come before the physical accuracy of the model. This philosophy is very close to the requirements for a light geological deformation algorithm presented in this paper.

A good review of the physically-based deformation algorithms is presented in [START_REF] Nealen | Physically based deformable models in computer graphics[END_REF]. They can be split into three main categories:

Continuum Mechanics Based Methods: these methods use different partitions of the space to solve the physical equations governing the deformation phenomenon described by continuum mechanics. The principal methods in this category are the Finite Difference Method (FDM), the Finite Element Method (FEM) and the Boundary Element Method (BEM). As these methods rely on actual physical laws, they generally converge towards the actual solution when the computation support is sufficiently fine.

Simplified Medium Methods: they approximate the deformed space with elements that deform with a simplified behaviour. For example, Mass-Spring Systems use a set of masses linked together by springs whose behaviour are described by very simple equations. Discrete Element Methods (DEM) also use this kind of simplification as they approximate the model by a set of particles interacting with each other by contact forces [START_REF] Burbidge | Numerical models of the evolution of accretionary wedges and fold-and-thrust belts using the distinct-element method[END_REF][START_REF] Scholtès | A DEM model for soft and hard rocks: Role of grain interlocking on strength[END_REF][START_REF] Vidal-Royo | The roles of complex mechanical stratigraphy and syn-kinematic sedimentation in fold development: insights from discrete-element modelling and application to the Pico del Águila anticline (External Sierras, Southern Pyrenees)[END_REF]. As these methods use an approximation of the physical laws, there is no guarantee they converge towards the physical solution.

Space Deformation Methods: these methods decouple the problem of computation of the deformation from the representation of the object to be deformed. The first space deformation methods were mainly driven by user manipulation [START_REF] Gain | A survey of spatial deformation from a usercentered perspective[END_REF], without constraints on the physical accuracy, e.g. Free Form Deformation (FFD) [START_REF] Barr | Global and local deformations of solid primitives[END_REF][START_REF] Sederberg | Free-form deformation of solid geometric models[END_REF] and Simple Constrained Deformation (SCODEF) [START_REF] Borrel | Simple constrained deformations for geometric modeling and interactive design[END_REF]. Subsequently, additional behaviours have been added to the space deformation approaches, e.g. by controlling the volume variations [START_REF] Hirota | Fast volume-preserving free-form deformation using multi-level optimization[END_REF][START_REF] Teschner | A versatile and robust model for geometrically complex deformable solids[END_REF][START_REF] Von Funck | Vector field based shape deformations[END_REF].

Two characteristics appear to be very important in the context of this paper: the space deformation, which meets the Space Consistency and Adapted Scale requirements, and the physical approximation, which meets the Interactivity over Physical Accuracy requirements. [START_REF] Botsch | Adaptive space deformations based on rigid cells[END_REF] present an algorithm which shows these characteristics together with very interesting robustness capabilities. As it is based on rigid elements and deforms the structures that are embedded into its elements, we refer to this algorithm as Reed (Rigid Elements Embedded Deformation). This acronym is a metaphor for the robustness of this model (it bends, but does not break).

Presentation of the deformation model Reed

Reed is a deformation algorithm based on the optimisation of the position of rigid elements with respect to a cost function E, similar to linear mechanical energy. The main steps of this algorithm are (Fig. 1):

1. Encapsulating the objects to be deformed into rigid elements (Fig. 1B). 2. Applying boundary conditions on some elements (Fig. 1C). 3. Iteratively optimising the rigid motion of the elements by minimising E (Fig. 1D). 4. Applying the interpolated displacements to the embedded objects (Fig. 1E). Most of the principles and equations presented in this section are adapted from [START_REF] Botsch | Adaptive space deformations based on rigid cells[END_REF]. However, the given details correspond to the developments of our implementation, whose formulation slightly differs from [START_REF] Botsch | Adaptive space deformations based on rigid cells[END_REF]. The principal difference relates to the description of the rigid element displacements. We consider the rotation components as a local rotation around the centre of the elements instead of a global rotation around the origin. This better separates translation and rotation increments and simplifies the process of computing rigid motion from the optimal linearised rigid displacements (see Section 3.3).

The encapsulation of deformed objects in rigid elements

The objects to be deformed are embedded in the rigid elements of the deformation model (Fig. 1B). All the displacement computations are made on these elements and the resulting interpolated displacement field is applied to the embedded objects at the end. This scheme simplifies the computation and guarantees the robustness of the approach for large deformations. The interpolation process is detailed in A.

When several objects are deformed, they are embedded in the same set of elements. Therefore, they are all deformed by the same displacement field, which makes it possible to preserve the relationships between the different deformed objects (provided that they are conformably discretised and that the discretisation is fine enough). The scale of Reed objects can be adapted to the scale of the deformation, independently of the scale of the geometrical details that are present in the structures to be deformed (Fig. 1), which reduces the number of elements required for the computation and allows fast computation.

The elements are displaced with a rigid motion (Fig. 2A), consisting of a rotation R around its centre c and a translation T:

M : R 3 -→ R 3 x -→ c + R(x -c) + T (1)
For convenience, all the elements used in this paper are cuboid, however Reed can also be applied to arbitrary polyhedral elements made to conform to geological features such as faults. The shape of the elements is only shown to represent the volume of initial space covered by each element for the energy definition (Section 3.2). They have no "space-filling" property.

Most of the computations performed on these elements require information about neighbour relationships between elements. The connectivity between the elements is given explicitly and is kept during the deformation process, which avoids the need to re-evaluate each element's neighbourhood. However, Reed elements do not have an explicit neighbourhood pattern. The number of neighbours per face is not fixed. This allows for locally refining the discretisation by splitting some elements, which is useful for locally adapting the elements to the embedded objects.

The cost function coupling the elements

The rigid elements are linked together by a cost function (Fig. 3). To stay close to the notations of [START_REF] Botsch | Adaptive space deformations based on rigid cells[END_REF], this cost function is referred to as the coupling energy E ij , even if it does not have the characteristics of energy in the strict sense.

E ij is defined as the integral, over the pair of elements C i , C j , of the difference between the displacements M i and M j , weighted by w ij [START_REF] Botsch | Adaptive space deformations based on rigid cells[END_REF]:

E ij = 1 2 w ij Ci∪Cj M i (x) -M j (x) 2 dx (2)
The weight w ij for two adjacent elements depends on their respective volumes, V i and V j , on the contact area S ij and on the distances h i and h j between each centre and S ij . It also allows us to introduce stiffness variations by weighting the contribution of certain elements or connections to the global energy. This is obtained by multiplying the geometrical factors by a dimensionless stiffness weighting w st :

w ij = w st ij S ij h i + h j V i + V j (3)
where h i and h j are the distances between the contact surface and the element's barycentre. The total energy E of the deformation model is obtained by summing all the coupling terms E ij .

The Reed energy formulation shares some interesting common points with the restoration algorithm based on rigid triangles proposed by [START_REF] Gratier | Restoration and balance of a folded and faulted surface by best-fitting of finite elements: principle and applications[END_REF]. The optimal position of the triangles are obtained by minimising the difference of displacement between the vertices of each triangle and those of the hole left by its neighbours. Reed's criterion slightly differs as it directly minimises the difference of rigid motion between each element and its neighbour.

The optimisation of the rigid element displacements

The deformation is driven by Dirichlet boundary conditions. These conditions are implemented by displacing elements and partly or fully fixing their displacement parameters. Similarly, as in [START_REF] Botsch | Adaptive space deformations based on rigid cells[END_REF], the displacement of the other elements are then optimised by a Newton-like iterative process to minimise the global energy E (Fig. 4). At each step k of this process, two kinds of displacements are considered for each element: the current rigid motion M k and a displacement increment A k . The descent direction is obtained by finding A k such that E is minimised.

To overcome the non-linearity of E (Eq. 2), a linearised form of the displacement increment A is used instead of an actual rigid motion. The incremented displacements of an element C i applied to a point x is thus approximated by A i M i (x) (Fig. 4), where A i is a combination of a translation vector v and a rotation vector ω, defining a rotation around the current centre of the element (i.e. after it has been affected by M ):

A k i : R 3 -→ R 3 x -→ x + v i + ω i × x -M i (c) (4)
where × refers to the cross product, when applied to vectors.

The translation vector v i and rotation vector ω i are the unknown of the descent direction estimation process. For each element, they are gathered into an unknown vector U i , which is of dimension 6×1. A global displacement unknown vector U is constituted by assembling the U i of all the elements. These approximated displacements are substituted to M in Eq. ( 2) yielding the linearised coupling energy E * ij :

E * ij = 1 2 w ij Ci∪Cj A i M i (x) -A j M j (x) 2 dx (5)
Applying this equation to all the elements yields a second order polynomial of the global displacement vector U (detailed in B):

E * = U T Q U + R T U + C (6)
where Q, R and C are the polynomial coefficients of E * . The optimal incremented displacement is then obtained by finding the minimum of E * , by cancelling its first order derivative with respect to U :

QU + 1 2 R = 0 (7)
Without appropriate boundary conditions, this linear system is not well posed as an infinite number of solutions exists. They correspond to all the rigid transformations of the initial position of the elements. A minimum of 6 degrees of freedom have to be constrained to produce a unique solution. More constraints can also be applied in order to impose a deformation to the Reed object. These constraints are set by fixing some degrees of freedom, i.e. setting some values of U to 0 depending on the boundary conditions to be applied. It allows us to remove the corresponding rows and columns of Q and R in Eq. ( 6). This yields the free terms Q f and R f , which are respectively of dimensions n f ×n f and n f ×1, where n f is the number of non-constrained degrees of freedom.

To compute the rigid motions of the next step M k+1 , the linearised displacement increment A k is re-transformed into a rigid motion by computing a rigid rotation increment from ω: the rotation axis is the unit vector parallel to ω with the same orientation and the rotation angle is given by arctan ω . M k is then updated by adding this rigid rotation increment to its rotation part and v to its translation part.

Testing specifications of Reed

This section inspects the capabilities of Reed to meet the specifications given in Section 2.2.

Apply boundary conditions

Update and fix the displacement of constrained elements

Compute displacement increments

Find A by solving eq. 7

Update current displacements

Robustness to extreme deformation

As a first application example, we deformed a model of the Annot Formation (SE France), adapted from [START_REF] Salles | Contrôles structuraux en 3 dimensions de la sédimentation turbiditique dans les chaînes plissées: exemple des Grès d'Annot (Sud-Est de la France[END_REF]. It consists of a series of structural horizons with several faults and an unconformity, all represented by triangulated surfaces. Reed is applied to these structures with an increasing deformation intensity (Fig. 5). The results illustrate(1) the robustness to extreme deformation, (2) the conservation of the relationships between the different structural elements. The complete optimisation runs in about 6 seconds on a personal computer (2×2.66 GHz, 4 Go RAM) with about 1200 rigid elements.

This example also illustrates the algorithm's ability to locally refine the discretisation (Fig. 5D) to ensure a better deformation optimisation and interpolation. This refinement can be automated, for example by refining a cell when it has accumulated too much energy [START_REF] Botsch | Adaptive space deformations based on rigid cells[END_REF].

In a second application we deformed a cross-section containing six layers having different stiffness w st and with one refined layer (Fig. 6). They are laterally shortened by fixing the position of the elements located on the left extremity and by translating the elements of the right extremity. In their final state, the layers have been highly shortened to a degree that elements overlap. This configuration could cause the elements of other deformation algorithms to degenerate, which would make it impossible to continue computing deformations. With Reed, the elements are rigid and therefore cannot degenerate, which allows for additional iterations, which eventually results in a more reasonable deformation.

Also, even if the resulting configuration of the elements looks degenerate, there are no self-intersections of the deformed horizons (Fig. 6B). In extreme cases, the displacements obtained can produce inconsistent deformation of the deformed objects. This is the case if two neighbour elements invert their position or if non-neighbour elements collide. However, such a behaviour is highly penalised by the deformation energy. Such degenerate cases should then appear only when the applied constraints are inconsistent relatively to each other or inconsistent with the mechanical capabilities of the deformed objects. For example, natural extreme deformation may induce faulting, but our implementation is presently unable to dynamically introduce faults during the deformation.

The configuration observed on the final state of the deformation may look degenerate since many elements are penetrating their neighbour's volume (Fig. 6A). However, this configuration is not actually degenerate in the sense that neighbour elements do not invert their positions and non-neighbour elements do not collide. Consequently, no intersections occur in the deformed layers (Fig. 6B).

Another demonstration of Reed's numerical robustness is shown by removing the constraints applied to the displaced elements (Fig. 6C). The rigid elements return rapidly to their initial position. This illustrates the capability of the system to continue computing even with extreme deformation and to rapidly adapt to boundary condition changes introduced by the user, which is crucial for its robustness and interactivity and for providing visual feedback to the user when boundary conditions change.

Kinematical and mechanical consistency

First-order kinematical and mechanical consistency is another crucial specification. It is investigated by comparing the resulting deformation with the results of a finite element method. It is important to note that these two methods are not based on the same equations, so they may not converge to the same result. FEM is used as a reference because it re- lies on mechanics and thus theoretically converges towards an acceptable physical solution. It is interesting to observe that very similar results are obtained with Reed on a simple case (Fig. 7).

The deformed object is a sandstone beam with dimensions 30×3×3m and a linear elastic model (with non-linear Green strain tensor). It is deformed with FEM and Reed having similar constraints, i.e. the left wall is fixed and the top of the right wall is constrained to go down 10m.

The finite element computations were performed with the Gocad plugin RestorationLab [START_REF] Durand-Riard | Balanced restoration of geological volumes with relaxed meshing constraints[END_REF][START_REF] Muron | Méthodes numériques 3D de restauration des structures géologiques faillées[END_REF]. The beam is discretised by a tetrahedral mesh containing 1146 nodes. The FEM calculations are performed directly on this mesh whereas Reed's computations are performed on a coarser mesh made of 288 cubic elements. Different numbers of elements have been tested to illustrate the evolution of the computation time when increasing the number of degrees of freedom 1). This table indicates that Reed is more than 10 times faster on this example even for the same number of degrees of freedom (both implementations being scalar). Also, Fig. 7 shows that a few elements are sufficient to obtain a good approximation of physical deformations.

The displacements obtained by the two methods are very similar and they are both consistent with the theoretical mechanical behaviour of beams. For example, Fig. 7B shows that the lower and upper third of the beam accommodate the most work during the deformation. The object local dilatations are also similar, even if Reed tends to produce dilatations that are slightly higher than FEM in this example. This can be explained by the fact that the volume variation is not explicitly penalised in Reed. Also, in the present form, Reed cannot reproduce certain physical behaviour such as the Poisson effect (e.g. thinning of a beam being stretched). Finally, the differences between the displacements are relatively low, except in the left part where it is of the same order of magnitude as the displacement, because the norms of the displacement and the differences are both nearly 0. 5 Building, editing and restoring geological models with Reed

This section illustrates three examples of situations where Reed can be successfully applied. We first present Reed as an editing tool for globally deforming an existing model. Then we demonstrate how it can be used to progressively build a complex synthetic model with several geological events. Finally, Reed is also adapted to the purpose of restoration.

Structural model editing

The editing of geological objects represents a first kind of application in which Reed can be used either for helping to generate an adapted initial geometry or to update geomodel geometry when new data becomes available or to reflect subsurface uncertainty. This can be helpful while constructing the model, for example to add small or large scale folds that are poorly constrained by data. A common approach to modelling structures is to optimise the geometry of triangulated surfaces with respect to data constraints [START_REF] Caumon | Surface-based 3D modeling of geological structures[END_REF][START_REF] Mallet | Geomodeling[END_REF]. To this end, having a good initial geometry is generally a key point, particularly when the structures are highly folded (e.g. Fig. 8 in [START_REF] Mallet | Discrete modeling for natural objects[END_REF]). Reed can help in building these initial surfaces by rapidly folding simple geometries making the construction of structural surfaces easier.

Another application consists of locally perturbing the geometry of a structural model in order to explore the uncertainties about the structures, for instance to reflect velocity uncertainty in time-to-depth conversions (Fig. 8). In this context, keeping consistent geometrical relationships between the different structural objects can be tedious. Furthermore, the editing has to be relatively fast and involve few parameters to keep the dimensionality of the problem relatively low. This can be done by restricting the deformation parameters to the displacement of a few Reed elements. For example, in Fig. 8, the elements on borders are fixed while a vertical displacement is specified for a central element. Such deformation thus involves only one parameter which is the vertical offset. The deformation in this example takes 0.6s, allowing us to rapidly explore the space of uncertainties, for example by sampling the possible different models while varying the values of the structural parameters. The topological relationship between the different layers are maintained (faults and unconformities).

Forward modelling

Reed can also be used to build forward models (Fig. 9) similar to what was presented by [START_REF] Jessell | Structural geophysics: integrated structural and geophysical modelling[END_REF]. Using a pseudo-mechanical deformation tool helps us to obtain more flexible deformations than the geometric operators of [START_REF] Jessell | Structural geophysics: integrated structural and geophysical modelling[END_REF]. Indeed, with a kinematic approach, the displacement field has to be specified for the whole object to be deformed whereas a dynamic approach allows us to specify only the displacement of some key elements, the complete displacement field being then a result of the deformation model.

We generated a model with two successive sequences: two sedimentation events leading to horizontal layers, each followed by a folding event. The resulting rather complex model, with its associated deformation history, was built within a few minutes.

Our approach is promising for modelling basins and deformed belts with complex polyphase tectonics, by separating the different episodes. To this aim, an inverse modelling approach could be used with such a tool to explore the space of possible deformation scenarios and to help modelling geological structures with respect to these possible scenarios as in [START_REF] Cardozo | Trishear modeling of fold bedding data along a topographic profile[END_REF]; [START_REF] Cardozo | Optimized trishear inverse modeling[END_REF].

Restoration of large deformations

Compared to most restoration methods, Reed focuses on robustness and interactivity rather than on simulation accuracy. It allows us to run approximate restorations very quickly, without necessarily having to generate a volumetric 3D mesh. This is useful for instance when only a few structural surfaces have been built and one wants to rapidly generate a petrophysical model in restored space. It can also help to characterise possible tectonic scenarios in the early stages of the modelling process.

This requires us to adapt Reed's boundary conditions to the purposes of restoration (Fig. 10), the classical boundary condition in restoration being the flattening of the uppermost stratigraphic horizon of the structural model. With Reed, deformed objects, including stratigraphic horizons, are not directly involved in the deformation computation. As a consequence, this condition has to be transferred to the rigid elements embedding the stratigraphic horizon. This is similar to the restoration of implicit surfaces [START_REF] Durand-Riard | Balanced restoration of geological volumes with relaxed meshing constraints[END_REF].

The successive tasks involved in the implementation of this constraint in Reed are as follows:

1. The elements which intersect the stratigraphic horizon are detected. 2. For each intersected element, the intersection is approximated by the mid-plane of the points contained by the element (Fig. 10B). 3. The z value of the center of each intersected element is computed by adding the signed-distance to the midplane to derived elevation of the target horizon. 4. The rotation of each intersected element is computed such that the direction going from the center of the element to its projection on the mid-plane is vertical. The boundary conditions computed in this way are progressively applied to the model while optimizing the position of the free elements until the surface is restored.

Our experience shows that the conditions in z and in rotation may conflict when the curvature of the surface is too large relative to the size of the elements. Indeed, some approximations are introduced in the computation of the boundary conditions because of the use of local mid-planes instead of the complete surface. The mid-planes represent a simplification of the surface as they cannot take into account its local curvature. This effect can be limited by locally adapting the size of the elements to the geometry of the stratigraphic horizon. In the future these conditions could also be relaxed by introducing them as soft constraints into the model rather than locking the corresponding degrees of freedom.

The restoration with Reed has been tested on the Han-sur-Lesse anticline, Belgium [START_REF] Bonniver | Étude hydogéologique et dimensionnement par modélisation du « Système-Traçage » du réseau karstique de Han-sur-Lesse (Massif de Boine -Belgique[END_REF][START_REF] Henrion | Integrated characterization and modeling of cave network: application to the karstic aquifer of Han-sur-Lesse (Belgium)[END_REF]. The faults contained in the model are not taken into account because they represent small displacements and occurred after the anticline formation. Only the lower horizons have been restored because the other ones have very closed limbs, which complicates the boundary condition computation as several surface parts located on opposite limbs may intersect the same rigid elements. Refining the elements would be a solution to avoid this problem. Fig. 11 shows the horizons in their deformed position and their restoration with three different Reed objects with an increasing number of elements. The horizon flattening error is color-coded on the restored surface and the histogram of this error is also presented. For each Reed object, 7 iterations have been performed to ensure the convergence of the restoration. Calculation times are indicated given the number of iterations and the time per iteration. The total times are 1.52s for 360 elements, 3.08s for 594 elements and 29.4s for 2700 elements. The decreasing restoration error as the number of elements increases illustrates the impact of a better intersection computation when setting the restoration boundary conditions. However, the computation time grows rapidly, which limits the benefits of refining the elements. Refining the elements only locally, around the restored surface, would be a solution to avoid this problem while minimizing the restoration error.

Conclusion and perspectives

In this paper we define original specifications for a light geological deformation tool, which represent the main contribution of this paper. This kind of deformation algorithm should help geologists to better integrate geological deformation histories at the early stages of the modelling process. Namely, they could help the geomodellers to interactively deform, forward model and restore the geological structures, as illustrated in Section 5. These specifications are oriented towards interactivity and robustness during interpretive tasks. They place the physical accuracy on a second level, which introduces some interesting perspectives for the range of deformation tools for geomodelling. We propose Reed as a potential answer to these specifications, but other algorithms may also be investigated.

The fact that Reed relies on rigid elements brings an efficiency to this algorithm and a robustness to large deformations. The deformed objects are embedded into the rigid elements, which means that the deformation is computed at the scale of the rigid elements and transferred afterwards to the embedded objects. This allows us to decouple largescale deformations and small-scale geometric details. Adaptive refinement is another strength of this algorithm which allows us to adapt the size of the elements to the details of the embedded objects or to the degree of deformation. These capabilities make this method applicable even in the case of extreme flexural deformation (Fig. 5 and Fig. 6).

The examples provided in this paper indicate that Reed matches the desired criteria in terms of interactivity, robustness and, to some degree, physical accuracy (Section 4.2). In the case of 3D restoration, our experiments suggest that this method is a good complement to existing kinematic methods and elastic geomechanical methods. Indeed, it introduces a new volumetric way to unfold layers affected by possibly large deformations.

However, some limitations remain in the use of Reed for geological deformations. The mechanical behaviour of the rigid elements, induced by the coupling energy minimisation, promotes a pure bending deformation. In natural stratigraphic layers, the orientation of the layering results in an anisotropic mechanical response [START_REF] Durand-Riard | Handling natural complexity in 3D geomechanical restoration, with application to the recent evolution of the outer fold-and-thrust belt, deepwater Niger Delta[END_REF]. Another prospect for future development of Reed is therefore to adapt the energy formulation to anisotropic behaviour.

A third limitation comes from the fixed connectivity of the elements. In large displacement contexts, faults may the same limits to ease comparison).

be required to let the relationship between elements evolve during the deformation. Taking into account faults in the deformation model is also important for future work. This last limitation is certainly more challenging as it would require us to adapt the geometry of the elements cut by a fault, to introduce the fault as an internal object in the rigid element description and to define a cost function that depends on the displacement of elements that are not initially neighbours.

A last limitation is that, for the sake of performance, the algorithm does not check for self-collision. However, in the context of very large deformations, for example with recumbent folds, one fold limb can come in contact with another. This kind of situation could be detected in order to avoid the two limbs intersecting.

Despite these limitations, the examples presented in this paper illustrate the capabilities of Reed and its utility for large deformation and 3D geomodel editing. The application to restoration in a large deformation context is promising and needs to be further examined. Other kinds of applications can also be explored, for example, building a structural model together with its associated deformation could help the determination of fracture location or stress evaluation for well planning in an industrial context.

Reed may not be the only method to fulfil the specifications we defined. However, it highlights the interest of using methods which balance space deformation and physically realistic behaviour. As suggested by [START_REF] Natali | Modeling terrains and subsurface geology[END_REF], Reed has provided a working example highlighting the application of Computer Graphics methods for interactively modelling geological structures.

A Displacement interpolation process

In [START_REF] Botsch | Adaptive space deformations based on rigid cells[END_REF], the embedded objects are deformed by interpolating the translation of the center of the elements, with inverse distance weighting or with a sum of Radial Basis Functions (RBF). Here, we propose to use the whole rigid motion of each neighbour element by taking the rotation of the elements into account.

Fig. 12 illustrates this process on a 2D example with a bilinear interpolation:

1. The point to be displaced x is located with respect to the element where it is embedded and its neighbours. This is described by the local coordinates u and v, as shown in Fig. 12A. 2. The contribution of each neighbour is computed by applying their final rigid motion M i to the point x. 3. The interpolated final position M (x) is computed as a weighted sum of the neighbours contributions:

M (x) = (1 -u)(1 -v) M 0 (x) + u (1 -v) M 1 (x) + (1 -u) v M 2 (x) + u v M 3 (x)
For the final rendering of the deformed object, higher degree interpolation (e.g. cubic interpolation or tri-harmonic RBFs, as proposed by [START_REF] Botsch | Adaptive space deformations based on rigid cells[END_REF]) can be used to ensure smoother results.

B Linear system computation details

This appendix details the computation of the polynomial coefficients Q, R and C, in the expression of the linearised energy E * (Eq. 6). These coefficients are respectively a matrix, a vector and a constant. They are computed by applying the linearised increment of displacement (Eq. 4) to the coupling energy equation (Eq. 2). This equation is derived in its matricial form by applying the following operations.

First, the linearised increment of displacement is expressed with respect to the unknown, namely the linear velocity v and the angular velocity ω:

A i M i (x) = M i (x) + v i + ω i × M i (x) -M i (c i ) = M i (x) + v i + ω i × R i (x -c i ) + T i + c i -R i (c i -c i ) -T i -c i = M i (x) + v i + ω i × R i (x -c i ) = M i (x) + P i (x) U i with P i (x) = 1 0 0 0 ρ i (x) 2 -ρ i (x) 1 0 1 0 -ρ i (x) 2 0 ρ i (x) 0 0 0 1 ρ i (x) 1 -ρ i (x) 0 0

4 3 12 6

 412 uses of deformation computation for structural modelling . . . . . . . . . . . . . . 3 2.2 Specifications for a light geological deformation algorithm . . . . . . . . . . . . . . . . . . 4 2.3 Overview of physically-based deformation algorithms . . . . . . . . . . . . . . . . . . . . . Presentation of the deformation model Reed 5 3.1 The encapsulation of deformed objects in rigid elements . . . . . . . . . . . . . . . . . . 5 3.2 The cost function coupling the elements . . . 5 3.3 The optimisation of the rigid element displacements . . . . . . . . . . . . . . . . . model editing . . . . . . . . . . . . 12 5.2 Forward modelling . . . . . . . . . . . . . . . 12 5.3 Restoration of large deformations . . . . . . . Conclusion and perspectives 13

Figure 1

 1 Figure 1 Rigid element embedding deformation approach. The approach is illustrated for the editing of the geometry of a fluvial channel. (A) The initial geometry shows a channel with a straight path. (B) The channel is embedded inside a Reed object made of 14 cubic rigid elements. (C) Some constraints (red) are applied to fix the position of certain elements. Only the translations are fixed here, the elements are still allowed to rotate. (D) The translation and rotation of the elements are optimised by globally minimising the energy of the Reed object (Fig. 4). (E) The displacement field is applied to the channel giving its final geometry. Note how the initial high-frequency details (e.g. the one pointed by the black arrow) of the geometry have been preserved.

Figure 2

 2 Figure 2 Rigid element displacement. (A) Two neighbour elements Ci and Cj (bold squares) are displaced (thin squares). Each element is moved by a rigid motion M , made of a rotation R around its center c and a translation T. (B) The displacement of both elements are applied to the point x yielding the final positions Mi(x) and Mj(x). The distance ∆M (x) between Mi(x) and Mj(x)represents the difference of displacement between Ci and Cj at point x. This quantity is the basis of the formulation of the coupling energy between neighbour elements (Eq. 2). The dashed square represents the position of Ci if it follows the displacement of Cj.

Figure 3

 3 Figure 3 The energy linking the rigid elements. This illustration corresponds to the Reed object used to deform the synthetic channel presented in Fig. 1. (A) The coupling energy rendered per element faces. (B) The coupling energy rendered per element.

Figure 4

 4 Figure4 Optimisation process. The optimisation process consists in optimising the elements' position (see Fig.2) by progressively applying small linear displacement increments A, made of a translation v and a rotation ω. During the optimisation process, the cost function (Eq. 5) penalizes the difference of incremented displacement ∆A M(x) . The whole optimisation process is presented in the flow diagram (right). The rigid body motion M of each element is updated with optimal displacement increment until a convergence criterion is reached.

EFigure 5

 5 Figure 5 Folding geological structures. (A) The initial structural model (Annot Formation, SE France) [Salles, 2010]. Each surface represents the interface between successive layers. (B and D) The final position of the rigid elements (constrained elements highlighted by red borders). (D) Certain elements have been refined in the region of the hinge of the fold. (C and E) The deformation is applied to the structural model while maintaining the relationships between the different structural elements, even with extreme deformation (D and E).

Figure 6

 6 Figure 6 Robustness during extreme deformation. (A) A cross-section is laterally compressed by displacing the upper 3 right border elements (red) to the left until an extreme deformation is obtained. Some layers have a lower stiffness (stippled line) or are more refined which produces different strain localisation between the layers. (B) Six horizon traces, originally flat and centred on each layer, are deformed by the Reed object in (A). (C) When removing certain constraints, the system returns to its original configuration, which illustrates its robustness.

Figure 7

 7 Figure 7 First order mechanical validity of Reed, comparison with FEM. (A) The displacements obtained by FEM (yellow solid with tetrahedral mesh) and those obtained with Reed (green transparent cubes). The embedding cubes are positioned over the left and right walls of the solid to make these walls correspond to the center of the elements. (B) The normalised energy of the rigid elements in their optimal position. (C and D) The dilatation of the object obtained by Reed and by FEM. (E) The difference between the two dilatations (positive values mean the dilatation is larger in Reed than in FEM and vice versa). (F) The logarithm of the ratio of the relative norm of the displacement difference highlights the relation between the difference of displacement and the displacement.

Figure 8

 8 Figure 8 Structural model editing. (A) A model of the Annot Formation (SE France), from Salles [2010] consisting of a series of structural horizons with faults and an unconformity. (B) The structural objects are embedded in rigid elements. The elements are refined to fit the structures. (C) The deformation is computed with boundary conditions fixing the position of the bottom elements on the left and right borders and fixing a given vertical displacement for a central element (yellow borders). (D) Result of the deformation applied onto the initial structural model. The red surface represents the initial position of the base horizon.

AFigure 9 Figure 10 Figure 11

 91011 Figure 9 Forward modelling with Reed. The deformation algorithm is used to represent the tectonics involved in a scenario with several successive phases. (A) Modelling of the deposition of a first stratigraphic sequence in onlap configuration. (B) A first deformation event is constrained by the boundary conditions (yellow) which fix the NW and SE corners while displacing the NE corner to the W. (C) The deformation is interpolated onto the first stratigraphic sequence and results in a non-cylindrical anticline localised in the north part of the model. (D) A second stratigraphic sequence is deposited. (E) A second deformation event is computed with this time the NE and SW corners that are displaced toward each other and NW-SE rows of elements that are constrained vertically to force the appearance of two parallel cylindrical anticlines. (F) This deformation is applied to both sequences which superimposes the two deformation events. (G) A SW-NE cross-section shows the final geometry. The cross-section location is shown on (F).

Table 1

 1 Comparison of the execution time of Reed and FEM. The time in seconds correspond to the time observed for Reed and FEM to converge for the bending of a beam (Fig.7). The FEM computation uses the Gocad plugin RestorationLab[START_REF] Durand-Riard | Balanced restoration of geological volumes with relaxed meshing constraints[END_REF];

	Method Number of Elements Number of degrees of freedom time (sec)
	FEM	1146	3438	5.3
	Reed	11	66	0.03
	Reed	288	1328	0.22
	Reed	672	4032	0.33
	Muron [2005]. Reed's result are obtained with our implementation in a Gocad plugin. All the times was obtained on a personal computer
	(2x2.66 GHz, 4 Go RAM).

Acknowledgements

The authors are grateful to the academic and industrial sponsors of the Gocad Consortium managed by ASGA (Association Scientifique pour la Géologie et ses Applications, www.gocad.org) for funding this work, and to Paradigm Geophysical for providing the Gocad software and API. This work was performed in the framework of the « Investissements d'avenir » Labex RESSOURCES21 (ANR-10-LABX-21).

and

where P i (x) is a 3×6 matrix and ρ i (x) a 3×1 vector.

The difference between the linearised displacement of two neighbour elements is expressed as:

The next step is to compute the squared norm of this term. The notations of the terms depending on x are simplified in the following equations to make it easier to read:

In this last equation, the term U T ij P T ij M ij is a scalar. As a consequence, it can be replaced by its transposed form, which yields:

E * is integrated over the volume of the two neighbours by a Gauss-Legendre integration scheme, which transforms the integral operator over each element into a sum over a given number of integration points. As the integrated term is a second order polynomial, it is integrated exactly with only 8 integration points for cubic elements. As compared to the analytical solution, this approach should yield exactly the same result and is easier to implement.

Finally, the integration process assembles the contribution of each pair of connected elements by gathering all the terms Q ij and R ij , which respectively forms Q and R, and yields Eq. ( 6).


Grimaldo Silva
email: jose.jgrimaldo@gmail.com

Anne-Hélène Olivier
email: anne-helene.olivier@inria.fr

Armel Crétual
email: armel.cretual@univ-rennes2.fr

Julien Pettré
email: julien.pettre@inria.fr

Thierry Fraichard
email: thierry.fraichard@inria.fr

Effective Human-Robot Collaboration in near symmetry collision scenarios

Recent works in the domain of Human-Robot Motion (HRM) attempted to plan collision avoidance behavior that accounts for cooperation between agents. Cooperative collision avoidance between humans and robots should be conducted under several factors such as speed, heading and also human attention and intention. Based on some of these factors, people decide their crossing order during collision avoidance. However, whenever situations arise in which the choice crossing order is not consistent for people, the robot is forced to account for the possibility that both agents will assume the same role i.e. a decision detrimental to collision avoidance. In our work we evaluate the boundary that separates the decision to avoid collision as first or last crosser. Approximating the uncertainty around this boundary allows our collision avoidance strategy to address this problem based on the insight that the robot should plan its collision avoidance motion in such a way that, even if agents, at first, incorrectly choose the same crossing order, they would be able to unambiguously perceive their crossing order on their following collision avoidance action.

I. INTRODUCTION

With the increased interest in applications that rely on robots sharing spaces with people, it is important that robots are able to integrate seamlessly into these spaces. A robot being able to determine and replicate the manner in which people avoid a collision with each other could allow people to act naturally around the robot [START_REF] Vassallo | How do walkers behave when crossing the way of a mobile robot that replicates human interaction rules?[END_REF], potentially reducing cognitive load for people in the environment [START_REF] Carton | Measuring the effectiveness of readability for mobile robot locomotion[END_REF].

There have been several research works on collision avoidance between people during locomotor task. More specifically, previous works highlighted the collaborative nature of this task and that the relative collision avoidance contribution between people depends on crossing order [START_REF] Olivier | Collision avoidance between two walkers: Role-dependent strategies[END_REF] and that crossing order is determined at the start of the collision avoidance interaction between people [START_REF] Knorr | Influence of person-and situation-specific characteristics on collision avoidance behavior in human locomotion[END_REF].

However, in certain situations, the manner in which these collision avoidance contributions should be performed is ambiguous -what we call symmetric collision avoidance scenarios. In these scenarios people are unable to determine their crossing order and are, by consequence, unable to guarantee effective collaboration in terms of complementary 1, where people may incorrectly choose to pass each other in opposite sides which negatively affects their efforts to avoid collision.

Interpersonal coordination during collision avoidance has been shown to be mainly determined by situational factors such as properties of the environment, relative positions, speed and heading of people rather than their height, gender or personality [START_REF] Knorr | Influence of person-and situation-specific characteristics on collision avoidance behavior in human locomotion[END_REF]. However, planning collision avoidance motions solely based on a prediction of crossing order is insufficient to accommodate scenarios of uncertainty as in the absolute worst case of future collision, that happens in symmetric collision scenarios, choosing the most common crossing order is no better than chance because people would choose to cross first or last with equal probability.

In this context, our work focuses on how a robot should move among people, what we call a Human-Robot Motion (HRM) problem. More specifically, we would like to replicate the collaborative nature of motion adaptations between people to avoid collisions during locomotor tasks, particularly whenever the exact manner in which motion adaptations should be performed is ambiguous.

Our objective is to find a collision avoidance motion that accounts for potential collaboration to avoid a collision in the best case while also minimizing the chances of crossing order remaining unclear in the worst case, where agents incorrectly assume the same crossing order. In order to accomplish the above-mentioned goal, we first present a novel measure called Minimum Predicted Distance with Goal Heading (MPDH), in order to analyze data obtained from experiments involving people in situations with different levels of crossing order ambiguity during a collision avoidance task while walking. This measure is used to evaluate whether the motion adaptations of each agent are having an overall positive or negative contribution to collision avoidance. We then establish a relation between the optic variable named derivative of the bearing angle and the uncertainty over crossing order. Using this relation a collision avoidance strategy that minimizes the potential negative impact of crossing order mismatch is then engendered.

II. BACKGROUND

People are not regular obstacles, they react to the motion of other agents. Accounting for this reaction allows a robot to minimize its negative impact in a social environment.

The Social Force Model (SFM) can reproduce the manner in which people move and react to other static and dynamic obstacles [START_REF] Helbing | Social force model for pedestrian dynamics[END_REF]. This model can be summarized as representing internal motivations of pedestrians through the use of attractive and repulsive forces to elements in the environment.

Using SFM, a robot was able to predict how standing people that are inadvertently blocking its way would move to allow the robot to pass [START_REF] Cosgun | Anticipatory robot path planning in human environments[END_REF]. Similarly, Ferrer and Sanfeliu relied on SFM to predict the reaction of several people to a given robot motion plan so that a trajectory that minimizes the combined reaction of people could be found. [START_REF] Ferrer | Proactive kinodynamic planning using the extended social force model and human motion prediction in urban environments[END_REF].

Even though SFM has been extensively used to simulate crowds, it was empirically demonstrated that distance-based repulsive forces have difficulty reproducing individual pedestrian behavior [START_REF] Steffen | A modification of the social force model by foresight[END_REF].

Another approach, named Reciprocal Velocity Objects (RVO) [START_REF] Van Den Berg | Reciprocal n-Body Collision Avoidance[END_REF] can also be used to replicate the behavior of crowds [START_REF] Bera | Realtime multilevel crowd tracking using reciprocal velocity obstacles[END_REF]. However, RVO was not modelled to naturally reproduce some behaviors of people during collision avoidance, for instance, people were shown to sometimes reverse their crossing order when avoiding collision [START_REF] Knorr | Influence of person-and situation-specific characteristics on collision avoidance behavior in human locomotion[END_REF] even though that entails longer time to avoid collision [START_REF] Vassallo | How do walkers behave when crossing the way of a mobile robot that replicates human interaction rules?[END_REF].

Recent works used joint trajectories to represent both usual and reversed crossing order of people during collision avoidance [START_REF] Kuderer | Teaching mobile robots to cooperatively navigate in populated environments[END_REF], [START_REF] Kretzschmar | Socially compliant mobile robot navigation via inverse reinforcement learning[END_REF]. This allowed estimation of possible choices of crossing order for each agent and its associated trajectory. In contrast, our work is focused on near symmetric situations in which crossing order decision is ambiguous which can cause agents to choose the same crossing order. Our collision avoidance approach intends to reduce the negative impact of this ambiguity in the collision avoidance.

III. CHARACTERIZING A FUTURE COLLISION

SCENARIO BETWEEN PEOPLE Before presenting our approach for collision avoidance between robots and people, it is necessary to first characterize a future collision scenario. Our focus is on dyadic collision avoidance situations where agents have goal directed trajectories. In order to reproduce behavior of people in these situations, we introduce concepts associated to crossing order determination and a metric of collision avoidance progress.

A. Minimum Predicted Distance

In order to evaluate the consequences of crossing order uncertainty, a way to numerically represent its impact in the collision avoidance progress is fundamental. In our work, this progress is measured between the start and end of the interaction between people, which we represent as the bounds of the time interval in which people deviate from their desired velocity in order to avoid future collision.

Collision avoidance progress can be measured in terms of a predictive variable named MPD [START_REF] Olivier | Minimal predicted distance: A common metric for collision avoidance during pairwise interactions between walkers[END_REF]. This variable is presented as the distance where people would meet in the future in case no motion adaptations are performed and is defined as

MPD(t) = ∞ min l=t (p r (t) + v r (t) • (l -t)) -(p p (t) + v p (t) • (l -t)) (1)
where p r (t) and p p (t) represent, respectively, the current position of an agent r and p at time t. In the same manner, v r and v p represent, respectively, the current velocity of agent r and p and • denotes the euclidean distance. Thus, whenever people are at a certain distance from each other and MPD is below a certain threshold, people adapt their motion in order avoid collision -increasing the value of MPD. On the other hand, in (near) symmetry situations agents might incorrectly choose the same crossing order which might not increase or even reduce MPD.

B. Derivative of bearing angle in collision avoidance

The bearing angle denotes the angle between the heading of an agent and another obstacle (dynamic or otherwise). In [START_REF] Cutting | How we avoid collisions with stationary and moving obstacles[END_REF], the derivative of bearing angle was found to correctly predict future collision and also to be a strong indicator of crossing order.

Given this context, let r and p represent a robot and a person respectively. From the local coordinate space of the robot, the bearing angle of r with respect to p is defined as α r,p (t) = atan2(y, x) where x = x p (t) -x r (t) and y = y p (t) -y r (t). The derivative of α r,p (t) is defined as

αr,p (t) = y x 2 + y 2 dx + x x 2 + y 2 dy , (2)
where x r (t) and y r (t) represent the x and y position of agent r at time t.

The relationship between the derivative of bearing angle and crossing order was explained in terms of the visual stimuli [START_REF] Cutting | How we avoid collisions with stationary and moving obstacles[END_REF]. The first crosser sees the last crosser diverge from the center of its field of view (bearing angle goes away from zero) as the time to collision diminishes, while the last crosser sees the first crosser converging towards the center of its field of view (bearing angle approaches zero). These roles were shown in [START_REF] Cutting | How we avoid collisions with stationary and moving obstacles[END_REF] and [START_REF] Olivier | Collision avoidance between two walkers: Role-dependent strategies[END_REF] to affect the manner in which people avoid collisions with each other and were shown to be strongly correlated to the derivative of the bearing angle at the start of collision avoidance.

In order to formalize the decision of an agent based on its crossing order, while also accounting for situations in which crossing order is undefined, such as in headon collision scenarios, the concept of homotopy class is used. This concept was defined in [START_REF] Kuderer | Teaching mobile robots to cooperatively navigate in populated environments[END_REF], [START_REF] Kuderer | Featurebased prediction of trajectories for socially compliant navigation[END_REF] as the side in which agents pass each other. In any given dyadic collision avoidance situation at least one agent has to decide whether they will pass the other on the left or on the right side. Each discrete decision about whether to pass left or right of someone represents a homotopy class decision. In [START_REF] Kuderer | Teaching mobile robots to cooperatively navigate in populated environments[END_REF], the homotopy class of each distinct trajectory pair was calculated using

Θ r,p = t αr,p (t)dt (3)
where αr,p (t) represents the derivative of the bearing angle between full trajectories of r and p. In this formulation, agents that pass each other on their right side will generate a Θ r,p = π while passing each other on the left side results in Θ r,p = -π. This formulation is elegant as its value is independent of the duration of collision avoidance. However, it is important to note that although Θ r,p correctly associates the choice of homotopy class of two actual (or predicted) trajectories, to complement these results [START_REF] Kuderer | Teaching mobile robots to cooperatively navigate in populated environments[END_REF] also built a probability distribution that is able to predict the most likely choice of homotopy class for any given person in the environment.

Our hypothesis is that predicting the probability of each distinct homotopy class decision is not sufficient. Understanding the underlying causes and consequences of this ambiguity allows one to explicitly accommodate the possibility of agents choosing actions that, when considered in conjunction with the robot's action, do not contribute to collision avoidance.

The concept of crossing order can then be used as a basis to predict, at any given time t, the homotopy class decision of an agent. However, attempting to directly utilize this information to determine collision avoidance motions is insufficient. The reasons for this can be contextualized given a direct representation of such predictor, presented as

Θr,p (t) =                +π or -π, if αr,p (t) = 0 -π, if α r,p (t) ≥ 0 and αr,p (t) > 0 -π, if α r,p (
t) < 0 and αr,p (t) > 0 +π, if α r,p (t) ≥ 0 and αr,p (t) < 0 +π, if α r,p (t) < 0 and αr,p (t) < 0 (4) which can then be simplified to

Θr,p (t) =      +π or -π, if αr,p (t) = 0 -π, if αr,p (t) > 0 +π, if αr,p (t) < 0 (5)
From this formulation, it is possible to verify that ambiguity could occur in two scenarios: whenever α r,p (t) ≈ 0 and αp r (t) ≈ 0, that is, a future (near) head-on collision, and in the (near) symmetric cases when α r,p (t) ≈ 0 and αr,p (t) ≈ 0. In both of these scenarios each person is met with a situation where +π and -π are possible solutions. As a consequence, in the space of solutions of these ambiguous scenarios, denoted as Θr,p (t) × Θp,r (t) = {{π, π}, { -π, -π}, {π, -π}, { -π, π}}, people may try to pass each other on the same side or may incorrectly try to pass each other in opposite sides, which will negatively affect their effort to avoid future collision e.g. decrease MPD. In both ambiguity scenarios, the uncertainty over crossing order depends on the value of α, our objective is to predict the uncertainty associated with this boundary, evaluate the potential consequences of this uncertainty in generated motions and, finally, minimize the risk of collision while allowing for effective collaboration.

IV. CROSSING ORDER UNCERTAINTY FROM PEOPLE DURING COLLISION AVOIDANCE

In this section we evaluate empirical data obtained from collision avoidance situations between people. Our focus is determining the chance that at least one will move in a way that changes their crossing order as predicted by α.

In our collision avoidance dataset, described in details in [START_REF] Silva | Human inspired effort distribution during collision avoidance in human-robot motion[END_REF], we evaluate dyadic interactions between people when avoiding collisions with each other. The experiments were composed of over 450 runs divided about equally into five distinct crossing angles: 30°, 60°, 90°, 120°and 150°. Several filters were used into the data to remove cases with too much noise resulting in 312 actual experiments for analysis.

The start and end of the interaction are denoted as, respectively, t i and t f . As such, our main focus is on the cases where the crossing order changed within the interval [t i , t f ] and whether these changes are associated to the derivative of bearing angle.

A. Behavior of people in scenarios of uncertainty

In order to directly visualize situations where crossing order is misjudged within our dataset, we modify MPD formulation slightly to obtain a monotonically increasing function whenever crossing order is respected, named Minimum Predicted Distance with Goal Heading (MPDH) and defined as

MPDH(t) = ∞ min l=t (p r (t) + v des r (t) • (l -t)) -(p p (t) + v des p (t) • (l -t)) (6)
where v des r is the agent r desired walking velocity at this time if there were no obstacles. This formulation of MPDH is In this sense, it represents the collision potential in terms of its impact on an agent's ability to head towards his goal. As such, whenever crossing order is respected its value increases monotonically during collision avoidance.

This formulation, in contrast with MPD, depends on having the complete trajectories, as in our case the goal of each agent and their desired speed (before collision avoidance) can be easily determined from their trajectories. This measure was used in the result analysis section (see Sec. VI-B) and to evaluate collision avoidance trajectories in our dataset.

As shown in Fig. 2, MPDH was used to determine the impact of crossing order changes in the progress of the collision avoidance.

B. Determining uncertainty boundaries

Our approach relies on predicting the behavior of agents when they are avoiding each other even in case of ambiguous crossing order. In Fig. 3 we depict some possible trajectories for reactive collision avoidance systems in case crossing order mismatches happen. To complement these results, in Sec. VI, we compare this expected result with our approach.

The derivative of the bearing angle α is a strong predictor of crossing order, as shown in [START_REF] Cutting | How we avoid collisions with stationary and moving obstacles[END_REF]. However, the role of each agent in a collision avoidance situation is not always clear. In this section, our focus is on first establishing a relationship between α and the certainty in crossing order determination using our dataset. Our objective is to find what is the chance that in the interval between [t i , t f ] the predicted crossing order changes, that is, the probability of

Θ(t i ) = Θ(t) for some t ∈ [t i , t f ].
The results, fitted to a sigmoid S(z) = 4 indicate that, at t i , if α approaches zero the likelihood of the crossing order changing within a given collision situation increases. As an extra measure of certainty, to guarantee this is not caused by reconstruction error in the position or heading, we manually evaluated the situations in which crossing order changed and verified that in most cases collision avoidance actions for agents, in terms of speed and heading changes or lack thereof, when observed in conjunction, were initially detrimental to collision avoidance. Fig. 3: Reactive collision avoidance system whenever agents incorrectly choose the same crossing order. Point in red is the collision in case of no change in trajectory. Collision point in orange happens if the agents plan collision avoidance only once. Collision in the blue point when both agents continuously choose to cross last. Whenever agents continuously believe they cross first they reach a stable forward motion without ever colliding with each other or reaching their goal.

C. Establishing crossing order in symmetric scenarios

In our work, we attempt to establish effective collision avoidance between agents in symmetric scenarios in up to n decisions, with agents taking two decisions per second. Effective collision avoidance means that both agents choose different crossing orders to avoid each other.

To that end, it is important to calculate the probability of correct collaboration D that would allow for Λ confidence in resolving crossing order after n consecutive attempts. In our work we define 95% as an acceptable Λ value. To find the value of D we equate the chance of n consecutive failures with 1 -Λ, as

1 -Λ = [1 -D] n (7)
which yields

D = 1 - n √ 1 -Λ (8)
This means that, at each time step, the agent should attempt to guarantee that the chance of resolving crossing order is at least D. As seen in Sec. IV-B, crossing order certainty can be calculated using S( αr,p (t)), as such, at this time step in case D ≤ S( αr,p (t)) the optimal collision avoidance action would generally suffice. However, our approach also anticipates situations with ambiguous crossing order by controlling the derivative of the bearing angle in an attempt to guarantee that, in the next time step, the chance of collaboration will be, in average, enough to achieve the desired confidence over n steps.

Consider that the robot r decides to avoid the collision with an agent p by changing its velocity to v new r , this decision is made based on the current velocity of the other agent, which we denote as v p (t). This motion 1 , in case of no change in the behavior of p, would yield αnew r,p . Similarly, agent p would use a similar logic to obtain a possibly distinct αnew p,r . Based on Eq. 2, the combined effect of these individual velocity changes on the instantaneous value of the derivative of the bearing angle can be calculated using αnew r,p + αnew p,r . Whenever agents incorrectly choose the same crossing order αnew r,p and αnew p,r will have distinct signs i.e. one will be positive and the other will be negative. Our objective is to calculate the probability P that αnew r,p + αnew p,r in the next decision step will be equal or larger than a given threshold κ. This would guarantee a certain confidence that even if agents at first incorrectly choose the same crossing order the crossing order would not remain ambiguous at the next decision step.

The desired value of P(z) where z = αr,p (t) at current time step can be calculated by equating the probability of two incorrect collaborations at the current and subsequent decision step to the current 1-S( αr,p (t)) and the subsequent one where αnew r,p + αnew p,r

was not larger than κ, with 1 Do note that although αr,p(t) = αp,r(t), the values of αnew r,p and αnew p,r are not necessarily equal as in this case we are evaluating the new velocity of one agent against the current velocity of the other agent (and vice-versa). Fig. 5: Relationship between the current αr,p (t), the desired number of decisions n to start effective collaboration and L.

P(z) = 1 - (1 -D) • (1 -D) 1 -S(z) (9)
Our final steps are then to determine the threshold value κ and to determine the manner in which one can guarantee | αnew r,p + αnew p,r |≥ κ with average probability P(z). To this end, let X be a random variable uniformly distributed over the interval [0, L]. From this random variable we select two points X 1 and X 2 . The distance between these points, denoted as Y = |X 1 -X 2 |, has an average value of M = E[Y ] = L 3 (see appendix for details) where E[•] is the expected value. Due to the locally linear shape of its function, we assume that a sample from Y that is m units away from the mean approximately respects the equality P(M ) = P(M -m)+P(M +m) 2 , this would mean that

P(E[Y ]) = E[P(Y )]
. This approximation allows us to establish that with M = κ we can finally determine the value of κ using the inverse of the sigmoid function (logit) presented in Sec. IV-B, denoted as S -1 (w) = b + 1 a log( w 1-w ), as a function of the desired probability of collaboration in the next time step with κ = S -1 (P(z)).

With these elements we can choose motions that respect the relation | αnew r,p + αnew p,r |≥ κ with average probability equal to P(z) by sampling αnew r,p from a uniform distribution with interval length L = 3κ.

The relationship between n, αr,p (t) and L are shown in Fig. 5. As n increases the size of L decreases until it it reaches zero which signifies that the agent is guessing crossing order based on just its current αr,p (t). This means that whenever αr,p (t) = 0 both agents would be choosing crossing order with odds no better than chance unless L > 0.

D. Discussion

Our approach to tackle near symmetry scenarios during collision avoidance is using communication through motion. The agent samples from a uniform distribution of size L the value αnew r,p in relation to the person's current perceived motion to have Λ confidence that effective collision avoidance will start in n time steps. Given our results, both n = 2 and n = 3 were found achievable for a robot. However, depending on the crossing order and time to collision n = 2 sometimes imposed a confidence smaller than Λ.

Another potential limitation is a consequence of a short time to collision, which would impose limits on the range of motions that can still avoid collision and by consequence limit the range of derivative of bearing angles that can be explored. Given this situation and considering that people generally need at least two steps for motion adaptation [START_REF] Patla | Understanding the roles of vision in the control of human locomotion[END_REF], we chose to focus our approach in situations where time to collision is above one second.

V. COLLISION AVOIDANCE UNDER NEAR-SYMMETRY

In the context of dyadic collision avoidance, our premise is that each individual agent will do its best to avoid future collision while preserving crossing order. The range of actions an agent can perform is limited by their maximum linear and angular velocity. Considering these constraints, each individual agent will attempt to minimize the change in their desired velocity v des r (t) that can avoid the future collision. Each agent, at first, will do the best they can to avoid a future collision, this can be represented as

v * r = arg min v∈F p r v -v des r (t) (10)
where F p r is the set of velocities for r where MPD of r with respect to p is larger than the threshold for collision and also respects crossing order i.e. attempts to guarantee that Θp r (t i ) = Θp r (t) for all t ∈ [t i , t f ]. However, as was shown in Sec. IV, crossing order can be misjudged which can have a detrimental impact on collision avoidance. As such, first we define an additional constraint on v based on the αnew r,p it would generate, represented as

αnew r,p = R( α- r,p , α+ r,p ) (11)
where R( αr,p , α+ r,p ) is a random value in the interval defined by αr,p and α+ r,p , which represent, respectively, the lower and upper bound in possible values of the derivative of bearing angle for r with respect to p that can avoid collision in a given crossing order. To generate collision avoidance velocities that respect crossing order we use an optimization approach named Sequential Least SQuares Programming (SLSQP) [START_REF] Kraft | Algorithm 733: TOMP -fortran modules for optimal control calculations[END_REF] as it can handle any combination of bounds, equality and inequality constraints.

To guarantee the desired confidence, it is necessary that L ≤ α+ r,p -αr,p . It is only necessary to select a subset of the interval α+ r,p -αr,p with L length. Examples of possible collision avoidance motions with a specific αnew r,p , from a discrete subset of values within an interval of size L, are shown in Fig. 6.

Given perceived safety and kinematic constraints, it is possible that L ≥ α+ r,p -αr,p , in these cases the agent would still attempt to resolve the collision with reduced confidence. However, when time to collision is below one timestep with L ≥ α+ r,p -αr,p the agent yields as it does not have the required confidence for human-like collaboration. These situations are more frequent whenever agents have a time to collision below one second.

In situations where crossing order has already a confidence value higher than Λ, which means that L = 0, it is possible to preserve crossing order even when one agent cannot contribute to avoid a collision (or its contribution would be insufficient). For instance, as shown in Fig. 1b, in many situations the first crosser is unable to accelerate further or change heading in a manner that preserves crossing order and they are thus unable to contribute to collision avoidance. Nonetheless, in scenarios with clear crossing order the first crosser would be able to rely on the last crosser to avoid collision even without any cooperation.

In any similar case, it is fundamental that crossing order is unambiguous as otherwise the last crosser could incorrectly Positive results, where the value of αr,p -α+ r,p is larger than L, are marked in bold. As crossing order is chosen randomly it can affect the length of αr,p -α+ r,p in similar scenarios (crossing last usually allows for larger interval). In the case with 120°crossing angle a smaller n was used. perceive himself as first crosser. This would ultimately result in collision if both agents are unable to contribute to collision avoidance as first crosser.

VI. SIMULATED EXPERIMENTS

Evaluating whether agents can avoid collision in near symmetry scenarios requires analysis of data from the agents' motion in situations with ambiguous crossing order. To that end, several simulated experiments analysing different aspects of our approach. are described in this section.

In our experiments, a holonomic motion model was chosen for the robot but the change in heading direction θ r between decision steps is bounded as -π 2 < θr < π 2 (rad/s) during collision avoidance to allow for more predictable motions. Similarly, the maximum speed of both agents was set as v max r = v max p = 1.7 (m/s) and initial speed as s initial r = s initial p = 1.2 (m/s). Finally, to avoid sudden turns during collision avoidance, which would negatively affect motion predictability for people, robot motions are smoothed using cubic splines.

A. Impact of linear and angular constraints on nearsymmetry mitigation

Although we establish the condition L ≤ [ αr,p , α+ r,p ] for Λ confidence in effective collision avoidance, it is important to evaluate when, in practice, the robot can obtain a reasonable length of the interval [ αr,p , α+ r,p ] when the constraints in both v max r and θr are taken into account. To that end, several random variations of collision avoidance scenarios were simulated over two crossing orders. The results, calculated considering Λ = 95% and showcased in Table I, indicate that, when time to collision is small, the length of the interval [ αr,p , α+ r,p ] also decreases. In these situations, collision avoidance requires sharper changes in speed and heading, this means that a smaller range of motions that are still able to avoid collision can be achieved. Fig. 7: Comparison of ideal scenarios where people always respect crossing order to our approach where ineffective collaboration is possible. Our collision avoidance approach mitigated the negative consequences of such event. Continuous lines indicate ambiguous crossing order while dashed gray lines indicate crossing order is no longer ambiguous.

In our analysis it was also shown that in smaller crossing angles three decisions can be necessary to achieve confidence larger than at least 90%. In contrast, collision scenarios with crossing angle of 120°allowed us to utilize up to n = 2 while maintaining Λ confidence, even though the situations for both crossing angles were generated with the same random process, this indicates that smaller changes in heading and speed generated larger changes in the derivative of the bearing angle when crossing angle was higher.

B. Collision avoidance trajectories under near symmetry

In order to assess whether generated collision avoidance motions properly mitigate the potential negative impact of ambiguous crossing order, a dozen collision avoidance scenarios where generated for each of the four tested crossing angles, these are: 0°, 30°, 60°and 90°.

In order to establish a comparison baseline with our approach, the collision avoidance of RVO [START_REF] Van Den Berg | Reciprocal n-Body Collision Avoidance[END_REF] was chosen due to its extensive use in the literature. As can be seen in Fig. 8: Several examples of MPDH where agents misjudge crossing order. In the upper plot, the MPDH increased less efficiently during the period where crossing order is misjudged, in the lower plot crossing order being misjudged negatively affected MPDH (its value decreases). Each color represents a distinct collision avoidance situation. Fig. 7, RVO is only able to reproduce collision avoidance motions where crossing order is respected. Our approach accounts for the possibility that agents misjudge crossing order in near scenarios.

The negative impact of ineffective collaboration in both agents ability to head towards their goal can be visualized with MPDH. In Fig. 8, an incorrect choice of crossing order can cause MPDH to increase less efficiently or even decrease i.e. negatively affecting their efforts to avoid collision.

In most cases, effective collaboration started after at most three decisions. In under 7% of the evaluated cases, the nearsymmetry situation was not solved and both agents stopped when time to collision was below one second.

VII. CONCLUSIONS

This work tackles dyadic collision avoidance situations between a robot and a person when the choice of crossing order for effective collaboration during collision avoidance is unclear. Approximating the uncertainty around this choice allowed our approach to choose a collision avoidance motion in such a way that, even if agents initially chose to cross each other on different sides, the robot and the person were generally able to perceive the correct homotopy class decision on their next decision step.

As future work, an interesting direction would be to run experiments in a virtual reality environment with a robot and a person controlling an agent. This approach avoids the risk of actual collision between robot and person. Moreover, assessing the impact of personal factors, such as cultural norms [START_REF] Chattaraj | Comparison of pedestrian fundamental diagram across cultures[END_REF], can be an interesting direction.

APPENDIX

In this section we provide a justification 2 of why given two random variables X 1 and X 2 that are uniformly distributed 2 Empirical validation code is available at https://github.com/ jgrimaldo/demonstrations over the interval [0, L] their distance Y = |X 1 -X 2 | will be in average E[Y ] = L 3 . Suppose that given X 1 and X 2 we also choose a third random variable X 3 , also uniformly sampled from an interval of length L. The value of X 3 will be between X 1 and X 2 , in average, 1/3 of the time as they are equiprobable. This means that the distance between X 1 and X 2 needs to be, in average, a third of the length of the interval and as such can be calculated simply as E[Y ] = L 3 .

Fig. 1 :

 1 Fig. 1: Robot's perspective in situations of near symmetry.

  (a) Reconstruction of human motion from data obtained by infrared cameras. (b) Trajectories in crossing order mismatch situation. First crosser in black. (c) Non-monotonically increasing MPDH.

Fig. 2 :

 2 Fig. 2: Human trajectory comparison when crossing order is misjudged.

1 1+e

 1 -a(z-b) with parameters a = 39.936914 and b = -0.000037, are shown at Fig.

Fig. 4 :

 4 Fig. 4: Evaluation of 312 collision avoidance scenarios, divided into five partitions. Each partition represents the probability that the predicted homotopy class does not change at any point in the interval [t i , t f ].

  (a) Head-on future collision with different initial speed. Crossing order is not defined.(b) Both assume first crosser role. Future collision with 90°crossing angle. Constrained by maximum velocity.(c) Both assume last crosser role. Future collision with 90°crossing angle. Less constrained by minimum velocity.

Fig. 6 :

 6 Fig. 6: Samples of collision avoidance motions with misjudged crossing order. Each color in a dashed line indicates a specific choice of derivative of bearing angle. Solid red and black lines indicate, respectively, initial motion and motion with αr,p .

*This work is partially supported by the Brazilian National Counsel of Technological and Scientific Development (CNPq) and the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) funded by the French program Investissement d'avenir.