N Passat

C Ronse Lsiit

J Baruthio

J.-P Armspach

AUTOMATIC PARAMETERIZATION OF GREY-LEVEL HIT-OR-MISS OPERATORS FOR BRAIN VESSEL SEGMENTATION

Reliable segmentation of 3D magnetic resonance angiography (MRA) is fundamental for planning and performing neurosurgical procedures, but also for detecting vascular pathologies. We propose here a method for brain vessel segmentation based on mathematical morphology tools. This method, devoted to phase-contrast MRA (PC-MRA) performs vessel segmentation by applying an adaptive set of grey-level hit-or-miss operators on each point of the MR data. High level anatomical knowledge modeled by a vascular atlas is used in order to adapt the parameters of these operators (number, size, and orientation) to the current position. The method has been performed on 30 PC-MRA cases composed of both phase and magnitude images. The results have been validated and compared to segmented data obtained by applying a region-growing algorithm on the same database. They tend to prove that the method is reliable for brain vessel detection and additionnally provides information on vessel size and orientation without requiring any post-processing step.

INTRODUCTION

Phase-contrast magnetic resonance angiography (PC-MRA) is a technique frequently used to provide 3D images of cerebral vascular structures. Indeed, the availability of precise information about brain vascular networks is fundamental for planning and performing neurosurgical procedures, but also for detecting aneurysms and stenoses. Many classical image processing tools have already been applied to the case of vessel segmentation. Nevertheless, nearly all the proposed algorithms use very little a priori knowledge and do not process multimodal data. In a previous paper [START_REF] Passat | Atlas-based method for segmentation of cerebral vascular trees from phasecontrast magnetic resonance angiography[END_REF] we proposed a first attempt to use anatomical knowledge to perform a segmentation on PC-MRA bimodal data. In this paper, we propose a new algorithm based on grey-level hit-ormiss transform using structuring elements designed accord-Thanks to the Équipe-Projet Multi-Laboratoires Imagerie et Robotique Médicale et Chirurgicale (EPML #9 CNRS-STIC) for funding.

ing to high level anatomical knowledge. This knowledge has been modeled by creating a vascular atlas from a set of segmented PC-MRA. This paper is organized as follows. Section 2 reviews some of the previous vessel segmentation techniques. Section 3 provides details on PC-MRA acquisitions. Sections 4 and 5 describe the anatomical and theoretical background necessary to understand and justify the algorithm proposed in Section 6. In Section 7, the method is tested by comparing results with those obtained from a region-growing algorithm. Discussion and projects are presented in Section 8.

RELATED WORK

The vessel segmentation methods can be classified according to the kind of strategy they use. The very first proposed methods were based on filtering [START_REF] Du | Vessel enhancement filtering in three-dimensional angiograms using long-range signal correlation[END_REF]. Algorithms based on mathematical morphology (hysteresis thresholding in [START_REF] Gerig | Symbolic description of 3-D structures applied to cerebral vessel tree obtained from MR angiography volume data[END_REF], grey-level erosions and dilations in [START_REF] Cline | Combined connectivity and a grey-level morphological filter in magnetic resonance coronary angiography[END_REF], or grey-scale skeletonization in [START_REF] Yim | Grayscale skeletonization of small vessels in magnetic resonance angiography[END_REF]) and region-growing [START_REF] Zahlten | Portal vein reconstruction based on topology[END_REF] have also been proposed.

Among all the existing strategies, very few have been designed to take advantage of multimodal data. A method, proposed in [START_REF] Sanderson | Simultaneous segmentation of MR and Xray angiograms for visualization of cerebral vascular anatomy[END_REF], for cerebral vascular structures visualization uses both 3D MRA and 2D X-ray images. An algorithm has also been proposed in [START_REF] Passat | Atlas-based method for segmentation of cerebral vascular trees from phasecontrast magnetic resonance angiography[END_REF] where angiographic and non angiographic data are involved in an atlas-based region-growing algorithm. Nevertheless, the simultaneous use of images from different modalities, as the use of a priori anatomical knowledge, is quite unusual. The algorithm presented here, based on grey-level hit-or-miss transform [START_REF] Ronse | A lattice-theoretical morphological view on template extraction in images[END_REF] uses bimodal PC-MRA data and integrates a priori knowledge modeled by a vascular atlas.

PC-MRA IMAGING

PC-MRA is a non invasive 3D MR technique [START_REF] Dumoulin | Three-dimensional phase contrast angiography[END_REF] enabling to generate two different images during a same acquisition. The first data, named magnitude image, contains anatomi- cal, but no vascular information, like a classical MRI. The second data, named phase image, is obtained by only conserving the signal of moving structures. Thus, the obtained image contains the flowing blood signal, plus noise and artefacts. Sagittal slices of such magnitude and phase images are illustrated in Fig. 1. The 30 MRA used in this paper have been performed on a 1 Tesla whole-body scanner (Gyroscan NT/INTERA 1.0 T from Philips, gradient slope 75 T/m/s). The flow encoding sequence called T1FFE/PCA uses a TR of 10 ms and a TE of 6.4 ms. The acquired images have dimensions varying between 256 2 ×150 and 256 2 ×180 voxels.

ANATOMICAL KNOWLEDGE

The vessel segmentation methods generally use very little a priori knowledge. Indeed, they are based on three main hypotheses: vessels correspond to the voxels of highest value in MRA; vessels have a tubular shape; vascular networks are organized in a tree structure. Nevertheless, brain vessels present other important properties. First, all the parts of the head do not contain the same quantity of vessels (this assumption takes into account the fact that MRA is unable to detect vessels having a diameter lower than 1 mm). Second, the size and orientation of brain veins and arteries depend on their position. Then, in each brain or head area, homogeneous properties concerning vascular density and vessel size and orientation can be observed. This allows to create vascular atlases for modeling this knowledge and to use them to guide segmentation tools. In the following we propose to use such an atlas in order to guide grey-level hit-or-miss operators.

GREY-LEVEL HIT-OR-MISS TRANSFORM

The hit-or-miss transform is a classical tool for extraction of templates from binary images. It uses a couple of struc-turing elements (A, B) where A (resp. B) has to fit the object (resp. the background). The binary hit-or-miss operation (⊗) can then be defined by:

X ⊗ (A, B) = {p ∈ ε | A p ⊆ X ∧ B p ⊆ X c },
where X is a binary object of a given Euclidian or digital space ε (in our case, ε ⊂ N 3 ).

Any increasing operator can be extended to grey-level images by considering them as stacks of binary images. However, the hit-or-miss transform, as defined previously is not an increasing operation. Another definition ( ) has then been proposed by Ronse in [START_REF] Ronse | A lattice-theoretical morphological view on template extraction in images[END_REF]:

X (A, B) = {p ∈ ε | A p ⊆ X ⊆ B p },
where A and B are structuring functions (A ≤ B). The greylevel hit-or-miss operation can then be defined by:

F (A, B) = {i (p,t) | A (p,t) ≤ F ≤ B (p,t) },
where F is a grey-level function and i (p,t) is the impulse function defined by:

i (p,t) = t if x = p -∞ if x p .
In [START_REF] Ronse | A lattice-theoretical morphological view on template extraction in images[END_REF], it has then been demonstrated that:

[F (A, B)](p) = (F A)(p) if (F A)(p) ≥ (F ⊕ B * )(p) -∞ otherwise ,
where B * is defined by B * (p) = -B(-p).

Choosing two structuring elements A e and B e , and two grey-levels a and b, with a ≥ b, we can define the grey-level structuring elements A and B by:

A(p) = a if p ∈ A e -∞ if p A e , and
B(p) = b if p ∈ B e +∞ if p B e ,
and then obtain:

F A = (F A e ) -a
[F (A, B)](p) =          (F A e )(p) -a if (F A e )(p) ≥ (F ⊕ Be ) + a -b -∞ otherwise .
Using this definition is equivalent to compare, at each point, the minimum intensity a min of all points within A e and the maximum intensity b max of all points within B e . If a min ≥ b max + ab, then the point belongs to the transform.

METHOD

The proposed method is based on the use of the grey-level hit-or-miss transform. This mathematical morphology operator has already been used by Naegel et al. in [START_REF] Naegel | Using grey-level hit-or-miss transform for segmenting the portal network of the liver[END_REF] for the detection of the portal vein entry in liver CT data. A fixedsized operator was then described to fit this precise anatomical structure. In this paper, we propose to use the grey-level hit-or-miss transform to carry out a whole segmentation of the cerebral arterial and venous trees.

Shape of the structuring elements

In order to fit the searched vascular structures, two kinds of structuring elements have been choosen. The element A, assumed to model the vessels, is a discrete sphere of radius r A , while the element B, modeling the neighboring background, is a set of 6 points sampled from a circle of radius r B , having the same center than A and oriented according to angles (θ B , φ B ). A set of such structuring elements is illustrated in Fig. 2 (left picture). The choice of spheres instead of ellipsoids for vessel modeling is justified by their higher robustness in case of vessel tortuosity. Moreover, the use of points sampled from a circle instead of a whole discrete circle allows a better detection of the background even in case of vessel junctions or bifurcations, with a lower computation time.

Use of the atlas

The immediate approach consisting in applying every couple of structuring elements in each point of the studied MRA data (D i ), to obtain the segmented image (D s ):

D s = {p ∈ D i | ∃(r A , r B , θ B , φ B , a, b), (D i A)(p) ≥ (D i ⊕ B) + a -b},
leads to a prohibitive computational complexity. Since this strategy can easily require more than 10 10 operator applications for one MRA data, it is necessary to reduce the algorithm complexity. This can be done by choosing, for each voxel, a subset of structuring elements sufficiently small to allow an efficient computation time and sufficiently large to find the vascular structures. This choice is done by taking into account anatomical knowledge modeled by a vascular atlas, the creation of which is fully detailed in [START_REF] Passat | Cerebral vascular atlas generation for anatomical knowledge modeling and segmentation purpose[END_REF]. This atlas can be defined by:

A : I → [0, 1] × P(R) × P([0, π[×[0, π[) x → (d, T, O) ,
where P(X) = {Y | Y ⊆ X}, d is the probability to find a vessel, T is the set of possible vessel thicknesses and O is the set of possible vessel orentations at the position x ∈ I. In the following, (d, T, O) will also be denoted by

(A d (x), A T (x), A O (x))
. This atlas provides a statistical estimation of vessels parameters for each voxel of a reference image. It can then be used to determine subsets of structuring elements adapted to every voxel of the image to process.

Algorithm

The method takes as input a PC-MRA of the whole head, composed of both phase (P) and magnitude (M) images (see Fig. 1), the vascular atlas (A), and the reference PC-MRA magnitude image (M A ) associated to the atlas. The first step is the superimposition of M A on M. This non rigid registration provides a 3D deformation field, used to obtain from A a new atlas A de f correctly fitting the processed MRA.

The second step is the choice of the set S of voxels that will be processed by the hit-or-miss operators. This step is important for reducing the computation complexity since the vascular structures generally represent less than 3% of the image volume. Removing from S all the voxels x such as A d de f (x) = 0 enables to avoid processing the background or areas presenting no vessels such as the cerebellum or the skull.

The third step consists in determining, for all x ∈ S , a subset of structuring elements and to apply them on the phase image P. For each x ∈ S , this subset is created by choosing structuring elements (A, B) such as r A ∈ A T de f (x), (θ B , φ B ) ∈ A O de f (x), and r B ∈]r A , 2.r A ]. Since the structuring elements are discrete objects, only discrete diameters and orientations (orientations according to discrete lines) are chosen. It has to be noticed that a fixed low value has been chosen for ab, enabling to detect vascular structures presenting a low contrast with the background. All the operators (A, B) are successively applied on P at the position x until one (or no one) of them matches a structure. If a couple (A, B) detects a vascular structure, the point x is added to the result image.

This method provides a skeleton image of the vascular tree. In order to obtain a volumic object, for each point x detected by a couple of structuring elements (A, B), a dilation by A is finally carried out. An example of such a segmented vascular tree is illustrated in Fig. 2 (right picture).

RESULTS AND DISCUSSION

The method has been implemented in C++, and run on a computer using a 3 GHz Pentium IV processor with 2 GB of memory. The average computation time is then 48 minutes to process one PC-MRA (38 minutes for the registration step and 10 minutes for the segmentation).

The 30 case database described in Section 3 has been used for validation. The results obtained with the proposed method have been quantitatively compared to those obtained with a previously described region-growing algorithm [START_REF] Passat | Atlas-based method for segmentation of cerebral vascular trees from phasecontrast magnetic resonance angiography[END_REF]. In average, 76% of vascular structures are correctly segmented by the proposed algorithm. The comparison also shows that the proposed method generates 23% of false negatives (vessels segmented by the region-growing method but missed by the proposed algorithm) but only 1% of false positives. This false positive ratio is quite satisfactory. Indeed, the method is essentially based on shape and size of structuring elements, but not on threshold values. Then it is not sensitive to non-tubular high-intensity artefacts which are generally the main cause of false positives. Moreover, by definition, the couple of structuring elements that enables to detect a vessel point immediately provides information on the vessel size and orientation, then avoiding a postprocessing quantification step. Nevertheless, this method still presents weaknesses concerning the false negative ratio. Indeed, it only uses spherical structuring elements to detect vessels. Then, the segmented vessels present tubular shape, even though their real shape is not tubular (many brain vessels have ellipsoidal or even triangular cross sections). Moreover, the smallest vessel are hardly segmented. For such structures, replacing the spherical structuring elements by line segments could provide better results.

CONCLUSION

This paper has presented a novel brain vessel segmentation method, based on grey-level hit-or-miss operators. The size and orientation of the operators are guided by anatomical knowledge modeled by a vascular atlas. This algorithm has been applied on a database of 30 PC-MRA, providing promising results. The method already provides a reliable vessel detection, presenting few false positive errors by comparison to a previously proposed region-growing algorithm. It also enables to directly obtain information on size and orientation of the segmented vessels. Its main weakness remains the unability to correctly reconstruct the smallest vessels or the vessels having a non circular cross-section. This problem is inherent to the current use of fixed-shaped spherical structuring elements. Then further work will now consist in adding the structuring elements shape as a supplementary parameter to be determined in each point of the image during the segmentation.

Fig. 1 .

 1 Fig. 1. PC-MRA sagittal slices. Left: magnitude image slice. Right: phase image slice. Since images are generated from the same acquisition, they are perfectly superimposed.

  , and F ⊕ B * = (F ⊕ Be )b, with B = {-b | b ∈ B}, and finally:

Fig. 2 .

 2 Fig. 2. Left: Part of the family of the structuring elements used to carry out the grey-level hit-or-miss segmentation. Right: Segmented vascular tree obtained with the proposed method. The grey-levels correspond to the radii of the structuring elements which permitted to segment the points.


Lucas Roche

Ludovic Anish Monachan

Ludovic Saint-Bauzel
email: ludovic.saint-bauzel@sorbonne-universite.fr

An experimental

come

An experimental study of parameters influencing physical Human-Robot negotiation in comanipulative tracking task.

Lucas Roche * , Anish Monachan * and Ludovic Saint-Bauzel * Abstract-The present paper investigates the relationship between human and virtual partners' behaviors during a comanipulative task requiring negotiation. More precisely, the study is focused on the influence of the Virtual Partner's (which was designed in a previous study) propensity to take the lead during conflicts over its human partner behavior. The experimental design includes both Human-Human and Human-Robot trials. The Human-Human trials serve as a baseline condition and help assessing the importance of the partner's nature (human vs robot) knowledge in the subjects behavior. The Human-Robot trials are used to observe the effect of different tuning of the Virtual Partner on the dominance relationship within Human-Robot dyads. The results of the study show that Human behavior is consistent across Human-Human trials, and is not influenced by the knowledge of their partner's nature. Moreover, a simple tuning of the negotiation thresholds of the Virtual Partner allows to significantly modify the dominance of the human partner, although this modification has no lasting effects.

I. INTRODUCTION

Since the early concept of cobots [START_REF] Peshkin | Cobot architecture[END_REF], a significant progress in control, conception and safety has brought natural physical human-robot interaction closer to reality. Robotic devices can now smoothly interact with their environment, and react to some amount of unknown parameters. Robots are increasingly used to work alongside humans and to cooperate with them for numerous tasks in a wide range of applications, from industry to health-care [START_REF] Santis | An atlas of physical human-robot interaction[END_REF] [START_REF] Goodrich | Human-robot interaction : A survey[END_REF].

Communication is a crucial point for efficient cooperation, both with human and robotic partners. While visual feedback is generally the main tool used in the design of modern robots, for its simplicity and versatility, there are situations where visual feedback becomes impractical. In situations inducing a high cognitive load [START_REF] Miller | The magical number seven, plus or minus two: some limits on our capacity for processing information[END_REF], such as surgery or piloting a plane, the visual channel can become saturated and lose its efficiency [START_REF] Van Veen | Tactile information presentation in the cockpit[END_REF]. On the other hand, haptic feedback is scarcely used in current designs, although humans process it efficiently, both alone and in combination with other type of information [START_REF] Jones | Kinesthetic sensing[END_REF] [START_REF] Ernst | Humans integrate visual and haptic information in a statistically optimal fashion[END_REF]. Haptic feedback also plays a key role in dyadic situations where people perform task together while in physical contact. Current research on physical Human-Human Interaction (pHHI) shows that haptic communication occurs naturally between humans [8] [9] [10] [11] [12], and allows to increase the performance of dyads in comanipulative tasks [START_REF] Groten | The role of haptic feedback for the integration of intentions in shared task execution[END_REF] [13] [14] [START_REF] Takagi | Haptic communication between humans is tuned by the hard or soft mechanics of interaction[END_REF] [16] [START_REF] Roche | High stiffness in teleoperated comanipulation: necessity or luxury ?[END_REF].

The ability for humans to communicate through the haptic channel has inspired multiple research studies aiming at reproducing these abilities in robots for a better physical Human-Robot Interaction (pHRI) [START_REF] Aydin | A new control architecture for physical human-robot interaction based on haptic communication[END_REF][19] [START_REF] Maeda | Probabilistic movement primitives for coordination of multiple human-robot collaborative tasks[END_REF][21] [START_REF] Kucukyilmaz | Intention recognition for dynamic role exchange in haptic collaboration[END_REF]. Providing such a capability to robots could help robot designer to increase the acceptability and the clarity of the robots behavior. It is a key issue encountered when a robot is brought alongside the human to cooperate. In our previous work [START_REF] Roche | Implementation of haptic communication in comanipulative tasks: A statistical state machine model[END_REF], a proof of concept that it is possible to have a virtual partner mimicking a human partner in a negotiation situation without hindering the performance of the dyad, nor changing the role dynamic of the subjects. A limitation of this proof of concept was that the parameters of the virtual partner were fixed to a set of values that was mainly due to the mechanical limitations of the device. A new device has then been developed in order to have a wider range of possibilities in torques and measures [START_REF] Roche | The semaphoro haptic interface: a realtime low-cost open-source implementation for dyadic teleoperation[END_REF]. The aim of this paper is to explore the main parameters that could influence the dominance of the virtual partner in the dyad negotiation. Based on our previous observations, we hypothesize that the Leader/Follower behavior of humans can be modeled by an intrinsic time/force threshold for negotiation, which varies for each person. An experiment in pHRI settings is thus conducted to test the hypothesis and investigate the influence of the Virtual Partner's thresholds on its dominance, and the behavior of its human partners.

II. MATERIALS AND METHOD

A. The SEMAPHORO Haptic Interface

In order to study pHRI while controlling the information that the subjects can communicate, the choice is made to use coupled haptic interfaces and a virtual environment to realize the experimental task. The SEMAPHORO1 haptic interface [START_REF] Roche | The semaphoro haptic interface: a realtime low-cost open-source implementation for dyadic teleoperation[END_REF] is designed for low-impedance, high precision tasks. It is designed specifically for the study of lightweight dyadic teleoperation. The teleoperation controller implemented in the interfaces allows to transmit forces and displacements from one interface to another, recreating the sensation of physical contact from a distance.

It is implementing a 4 channels teleoperation controller is presented in figure 1, where the control (Cm) is a PD controller (P gain : 110; D gain : 1) and the gain transmission is proportional gain (C2 = 1). T d represents the delay in the controller loop. It is equal to one loop period since the controller of both interfaces is run by the same process, reducing the communication delay to a minimum. The controller is able to provide a high precision as, it is able to measure the handle displacement with a discrete step 0.035 mm width. The bandwidth of the system is 10Hz with a -3db cut-off. The maximum force that can be transmitted at the handle is 23N. The total rotational inertia of the device is J h = 2.96x10 -4 kg.m 2 which represents an point-mass equivalent to 46g at the tip of the handle.

B. Experimental Task: Comanipulative Negotiating Task

The experiment consists in a co-manipulative tracking task that two subjects have to complete as a dyad, with either a human or a virtual partner.

The subjects are physically separated by an opaque screen to prevent any visual clue about the actions of their partners, and wear audio headphones playing white noise to prevent any auditory clue.

A path (white line over black background) is scrolling down on their monitor, at a speed of 35mm/s. The subjects use the haptic interfaces described previously to control the position of a massless virtual object, represented on their screen as a cursor (see Figure 2). The cursor is the same for both subjects, as they share control over a single common virtual object. The subjects are asked to keep the position of the cursor as close as possible to the scrolling path. To further incite each subject to cooperate, they are told that their goal is to maximize the common performance of the dyad. Feedback about the common performance is given by the color of the cursor, which changes based on the distance between the closest path and the cursor (see Figure 3):

• Green if |X cursor -X P ath | < 5mm • Yellow if 5mm < |X cursor -X P ath | < 10mm • Red if |X cursor -X P ath | > 10mm
The path is composed of a procedurally generated succession of curves, divided in two categories:

• The "BODY" category is composed of sinusoidal-like paths of random directions but fixed duration. The purpose of these parts is to keep the subjects focused on the task between each of the studied parts. • The "CHOICE" category is the aim of the experiment: at fixed intervals, the path splits into a fork, imposing a clear choice to be made concerning the direction that the subjects need to follow (see Figure 3). Considering that the subjects can neither see nor hear each other, the only way they can come to an agreement about the direction to choose is to use either the visual feedback from the monitor, or the haptic feedback from the handles.

While the path's structure is strictly the same for both subjects, each subject is encouraged to follow a highlighted trajectory. During the CHOICE parts, subjects receive some information about which side they have to choose [START_REF] Groten | The role of haptic feedback for the integration of intentions in shared task execution[END_REF]; this information can differ, creating situations of agreement or conflict, distributed in three cases. This is done by highlighting one of the two paths of the fork (see Figure 3):

• SAME: Both subjects have the same information, no conflict occurring. • OPPO: Opposite information given to each subject, inducing a conflicting situation.

• ONE: Only one subject has the information. This condition forces the subjects to be ready to take initiative in case they are the only one having information about the path to choose. It is designed to discourage subjects from keeping a passive strategy all along the trials.

The subjects are informed about these choices and the different decision types beforehand. The trials last 115 seconds, and include 2 SAME choices (both directions), 4 ONE choices (all possibilities), and 10 OPPO choices (equally distributed in the two possible possibilities) for a total of 16 choices.

C. Virtual Partner

An algorithm is designed to define the behavior of the virtual partner during the comanipulative task. The algorithm design is based on observations about the human behavior from previous experiments in pHHI for the experimental task. The virtual partner behaves as if it controlled a simulated handle to realize the task, creating displacements and forces perceived by its human partner during the task. The algorithm has access to information about the target trajectory, the position of its simulated handle, the position of the cursor on the monitor, the effort transmitted through its handle, and the position of its partner's handle through the position of the cursor and its own handle. The basic

STARTING POSITION STARTING POSITION LEADER MODE FOLLOWER MODE

VP takes initiative Human takes initiative

Fig. 4. Schematic functioning of the algorithm. The algorithm is designed to let the human lead the movement as a default. In the absence of human initiative, the virtual partner engage the movement toward its own target.

behavior of the Virtual Partner is explained hereafter, and more in-depth explanations can be found in [START_REF] Roche | Implementation of haptic communication in comanipulative tasks: A statistical state machine model[END_REF].

In the BODY parts, the algorithm follows the path. When confronted to a CHOICE, the algorithm generates a minimum-jerk trajectory [START_REF] Flash | The coordination of arm movements: an experimentally confirmed mathematical model[END_REF] from its current position to the target position, based on the choice it has to make. The minimum jerk is following this equation :

X(t) = X o + (X f -X o ) * (10τ 3 -15τ 4 + 6τ 5 ) (1)
where X o is the origin position, X f the final position and τ = (t -t start )/Duration. The duration of the motion is generated from a normally distributed variable (N (µ, σ))

based on the average and standard deviation of the human behavior data. Two situations are then possible: if the human takes initiative2 before the starting time of the virtual partner, it lets the human lead, entering "Follower Mode". The virtual partner generates a new trajectory to follow. This new trajectory is based on a minimum-jerk model starting at the current position of the virtual partner and ending at the new target. If the human partner did not initiate a motion before the beginning of the virtual partner's trajectory, the virtual partner takes the initiative, entering "Leader Mode", and starts its planned motion.

Once the virtual partner has started a motion in "Leader Mode", an ability to negotiate is provided in case the human wants to contest the choice. The algorithm measures the part of the interaction force between the partners which is directed toward a change of trajectory (negative if the virtual object is currently on the right, positive if the virtual object is on the left), if this interaction force exceed a fixed force threshold F threshold for a duration superior to a fixed time threshold T threshold , the virtual partner switches to "Follower Mode" and generates a new trajectory to follow the human. This change in trajectory can happen multiple times if the conditions are met. Figure 4 illustrate the statemachine algorithm.

D. Measures

The present experiment is focused on the changes in behavior within the dyads. The data from the experiment is analyzed from the CHOICE parts of OPPO type only, in a time window of 2 seconds centered around the path forking. In OPPO choices, the dyads members have contradictory suggested trajectories, and therefore must negotiate in order to reach a common agreement on the path to take (left or right). The Dominance measure is calculated as the percentage of choices "won" by a subject, that is the number of times the dyad ended up choosing his/her/its suggested trajectory, divided by the total number of choices analyzed. More complex ways of describing the dominance have been used [START_REF] Groten | Predictability of a human partner in a pursuit tracking task without haptic feedback[END_REF][26] but such sophistication was not deemed necessary for our analysis.

In Human-Human dyads, the subject with the highest mean dominance across the trials is designated as the Leader, while his/her partner is referred to as the Follower. In some case, data analysis for the Leader and Follower can be done separately to search for differences in Dominance patterns.

E. Statistical analysis

Statistical analysis of the results is performed using repeated measure Analysis Of VAriances(ANOVAs) [START_REF] Davis | Statistical methods for the analysis of repeated measurements[END_REF] when there are 3 or more values of the independent variables, and with two-tailed student t-tests [START_REF] Senn | The first t-test[END_REF] otherwise. The student t-tests are corrected using Bonferonni correction. Results in the next sections are given with the following form: ANOVAs (F-value, p-value), t-tests (Bonferonni corrected p-value, Cohenns d coefficient for size-effect). In both methods, pvalue represents the probability to be in the NULL Condition (the same population in our case) given a statistical model and the analysis used, p is compared to the significant level of 5%.

III. EXPERIMENTAL PROTOCOL A. Experimental Conditions

Three experimental conditions are tested during the experiment:

• Known Human Partner (KHP) In this condition, two human subjects collaborate on the task, and know that they are paired with a human. The teleoperation controller ensure that the positions of the interfaces, and thus the cursors are identifical for both subjects at all time. This additionally permits the transfer of haptic information within the dyad. • Hidden Human Partner (HHP) : In this condition, the subjects are told they perform the task with the Virtual Partner, but are actually paired with a human partner.

Outside the subject's beliefs, the rest of the parameters are identical to the KHP condition. • Known Virtual Partner (KVP) : In this condition, the subjects are paired with the virtual partner and are correctly told so. The condition is further divided in multiple subconditions depending on the thresholds values chosen for the virtual partner. • Hidden Virtual Partner : This condition was tested in the previous work [START_REF] Roche | Implementation of haptic communication in comanipulative tasks: A statistical state machine model[END_REF] and has shown no statistical difference with KVP. Consequently, it will not be included in this protocol.

B. Negotiation Thresholds: Force and Time

The virtual partner used in the experiments includes a possibility for negotiation in cases of conflict with the human. When the robot is leading the trajectory and the human disagrees on the direction chosen, he/she can contest this choice by applying an interaction force in the opposite direction. If the applied interaction force is superior to a fixed F threshold during a fixed duration T threshold , the virtual partner yields and concedes the lead to the human.

In the previous experiments [16][17], the values of F threshold and T threshold were chosen as the average interaction forces and duration observed in the experimental data for human dyads, it was equivalent to KVP1. These thresholds however influence the behavior of the virtual partner, and thus it is reasonable to think that they can influence the behavior of the human subjects, and the role dynamic within the dyad.

The experimental design includes multiple values of F threshold and T threshold . From a preliminary evaluation, thresholds values are chosen so that the difference between them is perceptible enough during the trials. Since all combinations of these thresholds must be tested in the protocol, the maximal number of threshold values must remain low, in order to keep the experiment duration acceptable. Three values of F threshold are chosen, and two for T threshold , for a total of 6 combinations of KVP conditions. The chosen threshold values and following KVP conditions are shown in Table I.

C. Hypothesis and Experimental Design

The first hypothesis that is tested is that for any of the thresholds considered, increasing its value will increase the virtual partner dominance compared to the human subjects.

The second hypothesis is that the behavior of the subjects in human-human condition can serve as a baseline in the experimental task. The design must thus include trials in KHP condition. Moreover, in order to verify that the KHP condition is a reliable baseline, multiple trials must be included in order to verify the stability of the subject's behavior in human-human interaction.

The third hypothesis tested is that the a priori on the partner's nature has no influence on the human behavior. This hypothesis was partly validated in a previous experiment [START_REF] Roche | Implementation of haptic communication in comanipulative tasks: A statistical state machine model[END_REF], where no difference were found between trials where the subjects were paired with the virtual partner and knew it (KVP), and trials where the subjects where paired with the virtual partner while thinking they were paired with a human (HVP -Hidden Virtual Partner). The reverse situation has however not yet been tested : Is there a difference between trials where the humans are together and know it (KHP), and trials where the humans think they are with the virtual partner but are actually together (HHP)?

The fourth hypothesis tested is that the Leader and Follower subjects in human dyads (KHP/HHP) behave differently when paired with the Virtual Partner (KVP). It is expected that different profiles of subjects lead them to behave more or less dominantly when paired together, and that these differences can also be seen in their behavior while paired with identical Virtual Partners.

The fifth hypothesis is that the threshold values of the virtual partner have an influence on the human behavior on the following trial(s). For example, a subject paired with an extremely dominant partner will tend to behave in a more following fashion than his/her base behavior in the next trials, even if the partner changes. In order to test this hypothesis, the experimental design must include transitions from all the different KVP conditions to a human-human condition. For practical reasons, all subjects won't provide data for all these transitions, and the data sets will be acquired over multiple dyads, at the cost of a reduced statistical test power. The transition will be done from KVP to HHP rather than KHP in order to prevent potential influence of the a priori on the partner's nature, if the previous hypothesis happens to be false.

To summarize, five hypothesis are tested during the experiment:

• H1 : Increasing threshold values leads to increasing virtual partner dominance. • H2 : The behavior of human in human-human dyads is consistent and time invariant. • H3 : The a priori of the partner nature has no influence on the human behavior. • H4 : Follower and Leader in KHP conditions behave differently in KVP condition. • H5 : Interaction with a partner of different dominance levels has lasting effects on the human behavior.

D. Experiment Design

The protocol includes 6 different combinations of thresholds values in the KVP conditions. In order to counterbalance the potential carryover effects, a balanced latin square design is adapted. The balanced latin square gives six different possible sequences of the experimental conditions, meaning the subjects will be grouped in 6 different categories. We decided to acquire two trials for each condition, and thus the latin square sequence was repeated twice for each subject.

The KHP condition (baseline) is tested at the start of the experiment, and at the end, which will allow to observe potential fatigue effects, or changes in behavior during the experiment. The HHP condition is also tested twice, just after the first KHP condition, which provide data for the KHP to HHP transition, and just before the last KHP condition, which provide data for the HHP to KHP condition.

Since the second HHP trials are done at the end of the KVP tests, and that the latin square design makes each sequence end with a different condition, all the six transitions from KVP to HHP are represented in the protocol, but by different subjects. Additional precautions are taken during data analysis to counterbalance the fact that not all subjects provide data for all conditions in this test. The final experimental protocol is illustrated in Figure 5, along with the tested hypothesis.

The subjects are distributed in six categories, which follow similar protocols, except for the order of the KVP conditions. Each category is tested with five subjects, for a total of 30 subjects. The 30 subjects were distributed as 15 dyads (8 male-male, 7 male-female). All participants were righthanded and had no prior knowledge about the study, experiments or experimental set-up. Data were collected from each dyad after each trial in the experiment.

IV. RESULTS

A. Consistency of Human-Human dyads behavior

No statistical difference is found between the first and last KHP condition (p = 0.491, d = 0.13). The results indicates that the human dyads behavior is consistent and time invariant, hence that the KHP condition is a reliable baseline during for experiment, consistently with hypothesis H2.

B. KHP vs HHP

The leader won 73% of the conflicting choices in KHP condition. The leader won 67.7% conflicting choices in the HHP condition. No significant difference is found in the dominance between KHP and HHP (p = 0.1, d = 0.4). This result reaffirms that prior knowledge of partner's nature has no influence in the dominance, which is consistent with the hypothesis H3.

C. Influence of the thresholds values

The dominance results for the six KVP conditions can be found in Table II. Two way additive ANOVA on the six KVP conditions shows a main effect force ((F(2, 492) = 22.8, p < 10 -4 )) and time ((F(2, 492) = 3.556, p = 0.04)) thresholds variations, with no interaction effect (F(5, 492) = 0.42, p = 0.66).

Post-hoc analysis is conducted on the effects of force and time threshold variations separately. Test data is thus aggregated for the different threshold values. The virtual partner won 25.1% of the choices with the 1N force threshold, 40.7% with 2.5 N, and 54.3% with 6.25 N. The virtual partner won 36.9% of the choices for a time threshold of 0.2s, and 43.2% for 0.5s. Results once again show significance for the influence of force (F(2, 492) = 22.9, p < 10 -4 ) and time (p = 0.02, d = 0.34) thresholds. These results also suggest a stronger influence of the force threshold compared to the time threshold.

The virtual partner dominance significantly increases with the thresholds values, which confirms the hypothesis H1.

D. Differences between Leader and Follower dominance

Subject data is separated for Leader and Follower, and tests are performed to observe the influence of the thresholds for Leader and Follower subjects. The different dominance means are exposed in Figure 6.

Statistical analysis shows a significant influence of the force thresholds for both the Leader (F(2,296) = 7.591, p < 10 -4 ) and the Follower (F(2, 296) = 5.248, p = 0.005). On the contrary, no influence of the time threshold emerges for the Leader (p = 0.14, d = 0.29) or Follower (p = 0.97, d = 0.006).

Overall, there is a significant difference in dominance between the Leader and Follower, both in KHP (p < 10 -4 , d = 0.255), and in KVP (p < 10 -4 , d = 0.92). The subjects that are Leaders in KHP are significantly more dominant when interacting with the Virtual Partner. This is consistent with hypothesis H4.

E. Lasting influence of the interaction with the VP

For this test, subject data is aggregated according to the type of the last KVP condition realized (1 to 6, see table I), with 5 subjects in each groups. A one-way between subjects ANOVA is then realized on the results in the following HHP condition, with the type of VP tuning as independent variable.

The results show no significant difference between the 6 groups (F(5, 24) = 1.016, p = 0.43), indicating that interaction with different types of Virtual Partners do not induce changes in the behavior of the subjects on subsequent trials. This hypothesis H5 is thus invalidated.

F. Discussion

The results of the experiment show that subject behavior is consistent across trials in human-human condition (KHP), moreover, a priori on the partner nature doesn't seem to affect the subject behavior (KHP vs HHP). Consequently, it can be assumed that humans exhibit an intrinsic and subject dependent tendency to dominance in negotiation situations in pHHI. It is possible that humans naturally develop and express different levels of confidence and leading behavior, and that this behavior is consistent throughout time. This result is consistent with research on Social Relation Model [START_REF] Kenny | A social relations variance partitioning of dyadic behavior[END_REF], and could be validated by conducting a similar experiment while pairing the subjects in multiple dyads combinations. This was done for example in [START_REF] Groten | Predictability of a human partner in a pursuit tracking task without haptic feedback[END_REF], where Groten et al. found that a majority of the variability in dominance in dyads was subject dependent. The fact that interaction with the Virtual Partner does not have lasting influence on the following HHP trials further confirms the consistency of human behavior in negotiation phase.

The experimental data also shows that in pHHI, a certain imbalance exists in the dyads, as one of the subjects is systematically more dominant than the other. Additionally, humans that tend to be Leader in pHHI condition stay more dominant in pHRI than those who tend to be followers. Our results are still showing an evolution of the dominance level with both kind of persons. The results in the KVP conditions show that changing the force and time thresholds of the Virtual Partner manipulates the subjects behavior in reaction.

V. DISCUSSION AND CONCLUSION

This paper shows that the force threshold mainly and the time threshold are able to influence the decision in a comanipulative tracking task during human-robot negotiation. This is interesting for two reasons. Firstly, it demonstrates a simple way to modulate the Virtual Partner behavior in "physical Human-Robot Interaction" situation. Previous research has shown that the ability to dynamically change the role allocation is a crucial point for efficient comanipulation [START_REF] Jarasse | Slaves no longer : review on role assignment for human-robot joint motor action[END_REF] [START_REF] Abbink | Haptic shared control: smoothly shifting control authority?[END_REF], and changing threshold values in our model is a way to attain this result. Secondly, if it is possible to manipulate subject behavior through changes in force and time thresholds during negotiation with the Virtual Partner, it could be possible to model that behavior using the same thresholds.

Modeling human behavior in "physical Human-Human Interaction" negotiation through a combination of time and force thresholds during interaction would allow advances in the comprehension of human behavior. It would also allow to facilitate the integration of "physical Human-Robot Interaction" protocols by enhancing the communicative ability of "Cobots". This enhancement requires to explore different kind of message that could be implemented in a Virtual Partner. For instance, the ability to choose between multiple paths (from near to far) on only one side of the path tracking or on two sides with asymmetrical paths positions. Combining different messages can be the first stones to pave the road of a kinesthetic communication. Such communication can help the human to understand naturally the robot intention/message during a physical cooperation. In order to improve our Virtual Partners, we will need to take into account the possible inter-individual variations of behavior in order to ensure the best cooperation and thus the best performances for Human-Robot comanipulation.

Fig. 1 .

 1 Fig. 1. Implemented 4 Channels Controller for both kinesthetic interfaces

Fig. 2 .

 2 Fig.2. Description of the experimental setup: Each participant uses a one dof haptic interface to share the control over a virtual object. Visual feedback about the position of the object is given on their respective monitors as a cursor.

Fig. 3 .

 3 Fig. 3. Illustration of the different decision types: SAME, ONE and OPPO. The data about the choices is recorded from a 2s timezone around the path's fork (in red on the leftward figure). Visual feedback about the dyad's performance is given through the color of the cursor (from left to right : green, yellow, red).

Fig. 5 .

 5 Fig. 5. Experimental protocol and hypothesis tested.

Fig. 6 .

 6 Fig. 6. Dominance results for Leader and Follower in the different threshold conditions

TABLE I VALUES

 I OF F threshold AND T threshold FOR THE DIFFERENT KVP CONDITIONS.

	T threshold	F threshold	1	2.5	6.25
	0.2		KVP1 KVP3 KVP5
	0.5		KVP2 KVP4 KVP6

TABLE II VIRTUAL

 II PARTNER DOMINANCE FOR THE DIFFERENT KVP CONDITIONS.

	T threshold	F threshold	1	2.5	6.25
	0.2		18.3% 33.8% 52.3%
	0.5		27.3% 42.9% 54.9%

System for human-robot physical manipulation evaluation (Système d' Évaluation de la MAnipulation Physique HOmme-RObot).

Taking the initiative is here defined as engaging a movement of the handle resulting in a displacement of superior to

35% of the distance between the starting position and the target.


Kacper Pluta
email: kacper.pluta@univ-paris-est.fr

Pascal Romon
email: pascal.romon@u-pem.fr

Yukiko Kenmochi
email: yukiko.kenmochi@esiee.fr

Nicolas Passat
email: nicolas.passat@univ-reims.fr

Honeycomb geometry: Rigid motions on the hexagonal grid

Euclidean rotations in R 2 are bijective and isometric maps, but they generally lose these properties when digitized in discrete spaces. In particular, the topological and geometric defects of digitized rigid motions on the square grid have been studied. This problem is related to the incompatibility between the square grid and rotations; in general, one has to accept either relatively high loss of information or non-exactness of the applied digitized rigid motion. Motivated by these facts, we study digitized rigid motions on the hexagonal grid. We establish a framework for studying digitized rigid motions in the hexagonal grid-previously proposed for the square grid and known as neighborhood motion maps. This allows us to study non-injective digitized rigid motions on the hexagonal grid and to compare the loss of information between digitized rigid motions defined on the two grids.

Introduction

Rigid motions on R 2 are fundamental transformations in 2D image processing. They are isometric and bijective. Nevertheless, digitized rigid motions are, in general, nonbijective. Despite many efforts during the last twenty years, the topological and geometric defects of digitized rigid motions on the square grid are still not fully understood. According to Nouvel and Rémila, this problem is not related to a limitation of arithmetic precision, but arises as a deep incompatibility between the square grid and rotations [START_REF] Nouvel | Configurations induced by discrete rotations: Periodicity and quasiperiodicity properties[END_REF].

Pioneering contributions to two-dimensional digital geometry in the square grid were made in 1970s (see [START_REF] Kong | Digital topology: Introduction and survey[END_REF][START_REF] Klette | Digital Geometry: Geometric Methods for Digital Picture Analysis[END_REF] for a survey). The main reason for taking this Cartesian approach was that objects defined on the square grid are easy to address in the computer memory, as well as image acquisition devices, whose sensors are generally organized on square grids. Moreover, the square grid has a direct extension to higher dimensions. However, such a Cartesian framework suffers from fundamental topological problems and therefore one has to choose different connectivity relations for objects of interest and their complements [START_REF] Kong | Digital topology: Introduction and survey[END_REF], or use a cell-complex framework [START_REF] Kovalevsky | Finite topology as applied to image analysis[END_REF], for example.

On the contrary, the regular hexagonal grid (called hereafter the hexagonal grid) suffers less from these problems since it possesses the following properties: equidistant neighbors -each hexagon has six equidistant neighbors; uniform connectivity -there is only one type of connectivity [START_REF] Middleton | Hexagonal Image Processing: A Practical Approach[END_REF][START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]. However, operations such as a sampling of continuous signals with the hexagonal grids are often computationally more expensive than ones defined on the square grid. Indeed, a digitization operator on the square grid can be defined by a rounding function, while there is no such counterpart on the hexagonal grid. Nevertheless, a practically useful digitization operator on a hexagonal grid has been proposed by Her [START_REF] Her | Geometric transformations on the hexagonal grid[END_REF].

Motivated by aforementioned issues, we study the differences between digitized rigid motions defined on both types of grids. To understand these differences, we establish a framework for studying digitized rigid motions on the hexagonal grid at a local scale. Our approach is an extension to the hexagonal grid of the former study of digitized rigid motions of the square grid, which is based on neighborhood motion maps [START_REF] Pluta | Bijective digitized rigid motions on subsets of the plane[END_REF]. This enables us to study the structure of groups induced by images of digitized rigid motions. It also allows us to focus on the issue of the information preservation under digitized rigid motions defined on the two grids. We provide a comparison of the information loss between the hexagonal and square grids. In addition, we present a complete list of neighborhood motion maps for the 6-neighborhood in Appendix.

2 Digitizing on the hexagonal grid

Hexagonal grid

A hexagonal grid H is a grid formed by a tessellation of R 2 by regular hexagons of side length 1/ √ 3, as illustrated in Figure 1. Let us consider the center points of these hexagons, i.e. the points of the underlying lattice Λ = Z 1 ⊕ Z 2 (see Figure 1(a)) where

1 = (1, 0) t and 2 = -1 2 , √ 3 2 t
. We can then regard H as the boundary of the Voronoi tessellation of R 2 associated with the lattice Λ.

Digitized rigid motions

Rigid motions on R 2 are isometric maps defined as

U : R 2 → R 2 x → Rx + t (1)
where t ∈ R 2 is a translation vector and R is a rotation matrix.

According to Equation (1), we generally have U(Λ) Λ. As a consequence, in order to define digitized rigid motions as maps from Λ to Λ, we commonly apply rigid motions on Λ as a part of R 2 and then combine the results with a digitization operator. To define such a digitization operator on Λ, let us first define a digitization cell of κ ∈ Λ where

C(κ) = x ∈ R 2 ∀β ∈ B, ( x -κ ≤ x -κ + β ) ∧ ( x -κ < x -κ -β ) , (1, 0) -1 2 , √ 3 2 (a
B = { 1 , 2 , 1 + 2 } ⊂ Λ. The digitization operator is then defined as a function D : R 2 → Λ such that ∀x ∈ R 2 , ∃!D(x) ∈ Λ and x ∈ C(D(x)
). We finally define digitized rigid motions as U = D • U |Λ . Note that this definition of the digitization operator is rather theoretical but not computationally relevant. For readers who are interested in implementing digitization operators in the hexagonal grid, we advise the method proposed by Her [START_REF] Her | Geometric transformations on the hexagonal grid[END_REF].

Due to the behavior of D, that maps R 2 onto Λ, digitized rigid motions are, most of the time, non-bijective. This leads us to define a notion of point status with respect to digitized rigid motions, as in the case of the square grid [START_REF] Ngo | Topology-preserving conditions for 2D digital images under rigid transformations[END_REF]. Definition 1. Let λ ∈ Λ. The set of preimages of λ ∈ Λ with respect to U is defined as S U (λ) = {κ ∈ Λ | U(κ) = λ}, and λ is referred to as an s-point where s = |S U (λ)| and is called the status of λ. The non-injective and non-surjective behaviors of digitized rigid motions result in the existence of 2-and 0-points, respectively (see Figure 1(b)).

3 Neighborhood motion maps and remainder range partitioning

Neighborhood motion maps

In R 2 , an intuitive way to define the neighborhood of a point x is to consider the set of points that lie within a ball of a given radius centered at x. This metric definition actually remains valid in Λ and it allows us to retrieve the classical notion of neighborhood based on adjacency relations. Definition 3 (Neighbourhood). The neighborhood of κ ∈ Λ (of a squared radius r ∈ R + ), denoted N r (κ), is defined as N r (κ) = κ + δ ∈ Λ | δ 2 2 ≤ r . Remark 4. N 1 corresponds to the 6-neighborhood which is the smallest, non-trivial and isotropic neighborhood of κ ∈ Λ.

As stated above, the non-injective and/or non-surjective behavior of a digitized rigid motion may result in the existence of 2-and/or 0-points. In other words, given a point κ ∈ Λ, the image of its neighborhood N r (κ) may be distributed in a non-homogeneous fashion within the neighborhood of the image N r (U(κ)), r ≤ r , of κ with respect to the digitized rigid motion U.

In order to track these local alterations of the neighborhood of κ ∈ Λ, we consider the notion of a neighborhood motion map on Λ defined as a set of vectors, each representing the rigid motion of a neighbor.

Definition 5 (Neighborhood motion map). The neighborhood motion map of κ ∈ Λ with respect to a digitized rigid motion U and r ∈ R + is the function defined as

G U r (κ) : N r (0) -→ N r (0) δ -→ U(κ + δ) -U(κ)
where r ≥ r.

In other words, G U r (κ) associates to each relative position of a point κ + δ in the neighborhood of κ, the relative position of the image U(κ + δ) in the neighborhood of U(κ). Remark 6. The maximal squared radius r of the new neighborhood U(N r (κ)) is slightly larger than r, due to digitization effect. In particular, we have r = 4 for r = 1.

For a better understanding of G U r (κ), we will consider a visual representation of the G U r (κ) functions as label maps. A first-reference-map L r associates a specific color label to each point δ of N r (0) for a given squared radius r (see Figure 2(a) for the map L 1 ). A second map L U r (κ)-associated to G U r (κ), i.e. to a point κ and a digitized rigid motion U-associates to each point σ of N r (0) the labels of all the points δ ∈ N r (0) such that U(κ + δ) -U(κ) = σ. Each σ may contain 0, 1 or 2 labels, due to the various possible statuses of points of Λ under digitized rigid motions (see Remark 2). Figures 2(b-c) provide some examples.

Note that a similar idea was previously proposed to study local alterations of the neighborhood N 1 under 2D digitized rotations [START_REF] Nouvel | Configurations induced by discrete rotations: Periodicity and quasiperiodicity properties[END_REF], and local alterations of N r under 2D digitized rigid motions [START_REF] Pluta | Bijective digitized rigid motions on subsets of the plane[END_REF] defined on the square grid.

Remainder range partitioning

Digitized rigid motions U = D • U |Λ are piecewise constant, which is a consequence of the nature of D. In other words, the neighborhood motion map G U r (κ) evolves noncontinuously according to the parameters of U that underlies U. Our purpose is now to express how G U r (κ) evolves.

U |N 1 (κ) is then non-injective
Let us consider a point κ + δ ∈ Λ in the neighborhood N r (κ) of κ. From Formula (1), we have

U(κ + δ) = Rδ + U(κ). (2)
We know that U(κ) lies in a digitization cell C(U(κ)) centered at U(κ), which implies that there exists a value ρ(κ) = U(κ) -U(κ) ∈ C(0).

Definition 7. The coordinates of ρ(κ), called the remainder of κ under U, are the fractional parts of the coordinates of U(κ) and ρ is called the remainder map under U.

As ρ(κ) ∈ C(0), this range C(0) is called the remainder range. Using ρ, we can rewrite Equation (2) as U(κ + δ) = Rδ + ρ(κ) + U(κ).

Without loss of generality, we can consider that U(κ) is the origin of a local coordinate frame of the image space, i.e. U(κ) ∈ C(0). In such local coordinate frame, the former equation rewrites as U(κ + δ) = Rδ + ρ(κ). Still under this assumption, studying the non-continuous evolution of the neighborhood motion map G U r (κ) is equivalent to studying the behavior of U(κ + δ) = D • U(κ + δ) for δ ∈ N r (0) and κ ∈ Λ, with respect to the rotation angle θ defining R and the translation embedded in ρ(κ) = (x, y) ∈ C(0), that deterministically depends on (x, y, θ). The discontinuities of U(κ + δ) occur when U(κ+δ) is on the boundary of a digitization cell. These critical cases related to U(κ+δ) can be observed via the relative positions of ρ(κ), which are formulated by the translated hexagonal grid H -Rδ, more precisely

C(0) ∩ (H -Rδ). Let us consider H -Rδ for all δ ∈ N r (0) in C(0), namely H = δ∈N r (0) (H -Rδ). Then C(0) ∩ H subdivides
the remainder range into regions, as illustrated in Figures 3 and4 for r = 1, called frames. Compared to the remainder range of the square grid [START_REF] Pluta | Bijective digitized rigid motions on subsets of the plane[END_REF], the geometry of the frames is relatively complex (see Figure 4). From the definition, we have the following proposition.

Proposition 8. For any λ, κ ∈ Λ, G U r (λ) = G U r (κ) iff ρ(λ)
and ρ(κ) are in the same frame.

Generating neighborhood motion maps for r = 1

From the above discussion, it is plain that a partition of the remainder range given by C(0) ∩ H depends on the rotation angle θ. In order to detect the set of all neighborhood motion maps for r = 1, i.e. equivalence classes of rigid motions U |N 1 (κ) , we need to consider critical angles that lead to topological changes of C(0) ∩ H. Indeed, from Proposition 8, we know that this is equivalent to computing all different frames in the remainder range. Such changes occur when at least one frame has a null area, i.e. when at least two parallel line segments which bound a frame have their intersection equal to a line segment. To illustrate this issue, let us consider the minimal distance among the distances between all pairs of the parallel line segments. Thanks to rotational symmetries by an angle of π 6 , and based on the above discussion, we can restrict, without loss of generality, the parameter space of (x, y, θ) to C(0) × 0, π 6 . Then, we found that there exist two critical angles between 0 and π 6 , denoted by α 1 and α 2 , with 0 < α 1 < α 2 < π 6 . Note that the angles 0 and π 6 are also critical, and denoted α 0 and α 3 , respectively. Figure 4 presents three different partitions of the remainder range for angles θ 1 ∈ (α 0 , α 1 ), θ 2 ∈ (α 1 , α 2 ) and θ 3 ∈ (α 2 , α 3 ), respectively.

Based on this knowledge on critical angles, we can observe the set of all distinct neighborhood motion maps M r = U∈U κ∈Λ G U r (κ) where U is the set of all rigid motions U defined by the restricted parameter space C(0) × 0, π 6 . The cardinality of M 1 is equal to 67. It should be noticed that κ∈Λ G U 1 (κ) is constant: 49 for any U, except for θ = α i , i ∈ {0, 1, 2, 3}. Indeed, we have κ∈Λ G U 1 (κ) = 1 for θ = α 0 , 43 for θ = α 1 , 37 for θ = α 2 and 30 for θ = α 3 . Such elements of the set M 1 are presented in Appendix.

(0, 0)

H -R(1, 0) H -R 1 2 , √ 3 2 H -R(-1, 0) H -R -1 2 , √ 3 2 H -R 1 2 , - √ 3 2 H -R -1 2 , - √ 3
2 Fig. 3. The remainder range C(0) intersected with the translated hexagonal grids H -Rδ, δ ∈ N 1 (0) for rotation angle θ = π 12 . Each hexagonal grid H -Rδ is colored with respect to each δ ∈ N 1 (0) in the reference label map L 1 (see Figure 2

Eisenstein rational rotations

In this section, we study the images of the remainder map ρ with respect to the parameters of the underlying rigid motion. In the case of the square grid, it is known that only for rotations with rational cosine and sine, i.e. rotations given by Pythagorean primitive triples, the remainder map does have a finite number of images; furthermore these images form a group [START_REF] Nouvel | Configurations induced by discrete rotations: Periodicity and quasiperiodicity properties[END_REF][START_REF] Pluta | Bijective digitized rigid motions on subsets of the plane[END_REF][START_REF] Nouvel | On colorations induced by discrete rotations[END_REF]. When, on the contrary, cosine or/and sine are irrational, the images form an infinite and dense set: any ball of non-zero radius in the remainder range intersects images of the remainder map. We will show that, in the hexagonal grid case, a similar result is obtained for rotations corresponding to the counterparts of the Pythagorean primitive triples in the hexagonal lattice. These are called Eisenstein primitive triples, namely triples (a, b, c) ∈ Z 3 such that 0 < a < c < b, gcd(a, b, c) = 1 and a 2ab + b 2 = c 2 [START_REF] Gilder | Integer-sided triangles with an angle of 60 •[END_REF][START_REF] Gordon | Properties of Eisenstein triples[END_REF]. We shall say that the rotation matrix R is Eisenstein rational if

R =        2a-b 2c - √ 3b 2c √ 3b 2c 2a-b 2c        = a c 1 0 0 1 + b c        -1 2 - √ 3 2 √ 3 2 1 2       
where (a, b, c) is an Eisenstein primitive triple. We must note that any rotation matrix can be written in this form, with a, b, c real numbers, not Eisenstein triples (and not even integers in general).

To begin, we focus on the rotation part of a rigid motion and define a group G = ZR 1 ⊕ ZR 2 ⊕ Z 1 ⊕ Z 2 and its translation G = G + t. We state the following result. Proposition 9. The group G is a rank two lattice if and only if the rotation matrix R is Eisenstein rational.

Proof. First, we note that the density properties of the underlying group is not affected by affine transformations, i.e. a lattice (resp. dense group) is transformed into another lattice (resp. dense group). Here we consider

X = [ 1 | 2 ] -1 so that {Xκ | κ ∈ Λ} = Z 2 .
Then we obtain

Ř = XRX -1 = X        2a-b 2c - √ 3b 2c √ 3b 2c 2a-b 2c        X -1 = a c -b c b c a-b c , (3)
and study Ǧ = ŘZ 1 0 ⊕ ŘZ 0 1 ⊕ Z 1 0 ⊕ Z 0 1 , instead of G. The generators of Ǧ are given by the columns of the rational matrix B = Ř I 2 where I 2 stands for the 2 × 2 identity matrix. As B is a rational, full row rank matrix, it can be brought to its Hermite normal form

H = [T | 0 2,2 ]
. The problem of computing the Hermite normal form H of the rational matrix B reduces to that of computing the Hermite normal form of an integer matrix: c ∈ Z is the least common multiple of all the denominators of B; compute the Hermite normal form H for the integer matrix cB; finally, the Hermite normal form H of B is obtained by c -1 H . The columns of H are the minimal generators of Ǧ. Since the rank of B is equal to 2, H gives a base (σ σ σ, φ φ φ), so that Ǧ = Zσ σ σ ⊕ Zφ φ φ. As H gives an integer base, c Ǧ is an integer lattice. Finally,

G = ZX -1 σ σ σ ⊕ ZX -1 φ φ φ.
Conversely, let us prove that G is dense if (a, b, c) is not an Eisenstein primitive triple (up to scaling). q, pq) also lies in Ǧ, has norm less than 3 e e e (hence less than ε) and is linearly independent from e e e (the matrix Ř has no eigenvectors).

Consequently, b b b 1 has rational coefficients, and we may take a, b, c integers with gcd equal to 1. Since cos θ = 2a-b c and sin θ = √ 3b 2c , we conclude that these form an Eisenstein primitive triple.

In the next section, we focus on a subgroup of G obtained from the intersection with the remainder range, i.e. Ḡ = G/Λ and its translation by t (modulo Λ). Proof. Up to the affine transformation of (3), we may consider the quotient group G = Ǧ/Z 2 and we note that for a primitive Eisenstein triple (a, b, c), any two integers in the triple are coprime [12, p. 12]. Let us first give a characterization of G. From the proof of Proposition 9, we know that any element x ∈ Ǧ is a rational vector of the form ( q 1 c , q 2 c ). By definition, {x} ∈ G iff there exist n, m ∈ Z such that {x} = {nb b b 1 + mb b b 2 }, i.e. there exist integers n, m, u, v such that

q 1 c + u = n a c -m b c , q 2 c + v = n b c + m a-b c
or, equivalently, q 1 + uc = anbm, q 2 + vc = bn + ambm .

A linear combination of both lines yields directly bq 1aq 2 = c(-cmbu + av), hence bq 1aq 2 ≡ 0 (mod c) is a necessary condition. It is also sufficient. Indeed, let us suppose that bq 1aq 2 = kc, k ∈ Z. Then, since gcd(a, b) = 1, we know that the solutions to this Diophantine equation are of the form (q 1 , q 2 ) = (a, b) + kc(β, -α), where αa + βb = 1 (Bézout identity) and ∈ Z. Consequently, (

q 1 c , q 2 c ) = b b b 1 + (kβ, -kα) lies in Ǧ.
Moreover, we deduce that G is cyclic with generator b b b 1 . Finally, { b b b 1 } = (0, 0) implies a = uc and b = vc for some integers u, v. Applying the Gauss' lemma to coprimes a, c, we see that needs to be a multiple of c. Therefore | G| = c.

5 Non-injective digitized rigid motions

Non-injective frames of the remainder range

In Remark 2, we observed that some frames of the remainder range correspond to neighborhood motion maps which exhibit non-injectivity of the corresponding digitized rigid motion. The set of all the neighborhood motion maps M 1 , presented in Appendix, allows us to identify such non-injective zones i.e., unions of frames of the remainder range. For instance, the frames related to the neighborhood motion maps of the axial coordinates; (-2, 4), (-1, 4) and (0, 3) (see Figure 6 (top)), constitute a case of such a non-injective zone for rotation angles in (α 0 , α 1 ). Based on this observation, we can characterize the non-injectivity of a digitized rigid motion by the presence of ρ(κ) in these specific zones (illustrated in Figure 4(a)).

Conjecture 11. Let C 6 stand for the 6-fold discrete rotation symmetry group. Given U ∈ U and κ ∈ Λ, U(κ), has two preimages κ and κ + δ where δ ∈ N 1 (κ) iff ρ(κ) is in one of the zones c k (F ), c k ∈ C 6 , where F is the parallelogram region whose ver-

tices are: cos θ -1, cos θ √ 3 , 0, 1 √ 3 , 1 2 2 -cos θ - √ 3 sin θ , 1 6 √ 3 cos θ + 3 sin θ , 1 2 cos θ - √ 3 sin θ , 1 6 3 sin(θ) + 3 √ 3 cos(θ) -2 √ 3 .

Non-injective digitized rigid motions in square and hexagonal grids

In this section, we compare the loss of information induced by rigid motions on the square and the hexagonal grids. Indeed, we aim to determine on which type of grid digitized rigid motions preserve more information.

In accordance with the discussion in Section 4 and the similar discussion for the square grid in [START_REF] Nouvel | Configurations induced by discrete rotations: Periodicity and quasiperiodicity properties[END_REF][START_REF] Berthé | Discrete rotations and symbolic dynamics[END_REF], the density of images of the remainder map in the non-injective zones is related to the cardinality and the structure of Ḡ. On the one hand, when Ḡ is dense, it is considered as the ratio between the area of the union of all the non-injective zones k=1,...,6 c k (F ), and the area of the remainder range [START_REF] Berthé | Discrete rotations and symbolic dynamics[END_REF]. On the other hand, when Comparison of the loss of information induced by digitized rigid motions on the hexagonal and the square grids. The red and blue curves correspond to the ratios between the areas of noninjective zones and the remainder range for the square and the hexagonal grids, respectively. The "x" (resp. "o") markers correspond to the values of the information loss rate 1

-|U(S Z )| |S Z | (resp. 1 -|U(S L )| |S L | ) Ḡ forms a lattice, it is estimated as | Ḡ∩( c k (F ))| c , k = 1, . . . , 6
, where c is an element of a primitive Eisenstein (resp. Pythagorean) triple, i.e. | Ḡ| [START_REF] Berthé | Discrete rotations and symbolic dynamics[END_REF].

Nevertheless, to facilitate our study, we will consider the former ratio between the areas as an approximation of the information loss measure. Indeed, the rotations induced by primitive Eisenstein (resp. Pythagorean) triples are dense; one can always find a relatively near rotation angle such that c is relatively high. The area-ratio density measure can be then seen as a limit for the cardinality based ratio. Figure 5 presents the analytical curves of the area ratios for the hexagonal and the square grids with respect to rotation angles.

In order to validate this approximation, we also measure the loss of points for different sampled rotations and the following finite sets S L = Λ ∩ [-100, 100] 2 (resp. S Z = Z 2 ∩ [-100, 100] 2 ) in the hexagonal (resp. square) grid. In this setting, we use 1

-|U(S L )| |S L | (resp. 1 -|U(S Z )| |S Z |
) as the measure. They are plotted as well in Figure 5. We note that the experimental results follow the area-ratio measure that provides consequently a good approximation 5 . The obtained results allow us to conclude that digitized rigid motions on the hexagonal grid preserve more information than their counterparts defined on the square grid.

Conclusion

In this article, we have extended the framework of neighborhood motion maps to rigid motions defined on the hexagonal grid-previously proposed for digitized rigid motions defined on Z 2 [START_REF] Pluta | Bijective digitized rigid motions on subsets of the plane[END_REF][START_REF] Nouvel | On colorations induced by discrete rotations[END_REF]. Then, we studied the density of images of the remainder map and the structure of the induced groups. Finally, we have shown that the loss of information induced by digitization of rigid motions on the hexagonal grid is relatively low, in comparison with those on the square grid.

Our main perspective is to use the proposed framework for studying bijective digitized rigid motions and geometric and topological alterations induced by digitized rigid motions on the hexagonal grid.

Following a paradigm of reproducible research, the source code of the tool designed for studying neighborhood motion maps on the hexagonal grids is provided at the following URL: https://github.com/copyme/NeighborhoodMotionMapsTools Appendix: Neighborhood motion maps for G U 1 and their graph

In Section 3.3, we observed that there exist three topologically different partitions of the remainder range. They depend on rotation angles: θ 1 ∈ (α 1 , α 2 ), θ 2 ∈ (α 2 , α 3 ) and θ 3 ∈ (α 2 , α 3 ), as illustrated in Figure 4. From Proposition 8, we also know that each frame corresponds to a different neighborhood motion map. Figure 6 illustrates all the neighborhood motion maps of M 1 , with the different rotation angles: θ 1 ∈ (α 1 , α 2 ), θ 2 ∈ (α 2 , α 3 ) and θ 3 ∈ (α 2 , α 3 ).

The dual graph of the remainder range partitioning, in Figure 4, is considered. In this graph G 1 = {M 1 , E}, each node is represented by a neighborhood motion map (or a frame), while each edge between two nodes corresponds to a line segment shared by frames in the remainder range. Moreover, each edge is labeled with the color of the corresponding line segment (see Figure 3). We note that an edge color denotes the color of a point in the transition between the two corresponding neighboring motion maps. For instance, if there exists a red horizontal edge between two nodes, then we observe a transition of the red neighbor (δ = (1, 0)) between two neighborhood motion maps connected by this edge. We invite the reader to consider Figure 6 (top) and neighborhood motion maps of the indices (-1, 2) and (0, 2).

Note that, neighborhood motion maps in Figure 6 are arranged with respect to the hexagonal lattice and each can be identified thanks to its axial coordinates [START_REF] Middleton | Hexagonal Image Processing: A Practical Approach[END_REF]. We also observe that, in such an arrangement, neighborhood motion maps are symmetric with respect to the origin-the frame of index (0, 0). For example, the neighborhood motion map of index (-4, 3) is symmetric to that of the index (4, -3) (see Figure 6). r (κ) for any angle in (α 0 , α 1 ) (top), (α 1 , α 2 ) (bottom, left), and (α 2 , α 3 ) (bottom, right), visualized by the label map L U 1 (see Figure 2). Neighborhood motions maps are considered as graph vertices and linked by edges with respect to the adjacency relations of respective frames, i.e. each edge in the graph represents a side shared by two frames in the remainder range. The neighborhood motion maps which correspond to noninjective zones are surrounded by pink ellipses. The elements which are changed with respect to the rotation angles in the bottom figures are surrounded by black squares, while those which are not changed are faded. A high resolution version of the figure can be found via the URL: https://doi.org/10.5281/zenodo.820911

Fig. 1 .

 1 Fig. 1. A visualization of a part of the hexagonal grid (a) . The arrows represent the base of the underlying lattice Λ with the white dots representing its elements. (b) Examples of three different point statuses induced by a rigid motion U on the hexagonal grid: digitization cells corresponding to 0-and 2-points are marked by green lined , red dotted patterns and labeled with their status numbers, respectively. The white dots indicate the positions of the images of the points of the initial set Λ under U

Remark 2 .

 2 For any λ ∈ Λ, |S U (λ)| ∈ {0, 1, 2}. From the digitization cell geometry, we have that |S U (λ)| = 2 only when two preimages κ, σ ∈ S U (λ) satisfy κσ 2 = 1.

Fig. 2 .

 2 Fig. 2. (a) The reference label map L 1 . (b-c) Examples of label maps L U 1 (κ). (b) Each point contains at most one label: the rigid motion U |N 1 (κ) is then locally injective. (c) One point contains two labels: U |N 1 (κ) is then non-injective

Fig. 4 .

 4 Fig.3. The remainder range C(0) intersected with the translated hexagonal grids H -Rδ, δ ∈ N 1 (0) for rotation angle θ = π 12 . Each hexagonal grid H -Rδ is colored with respect to each δ ∈ N 1 (0) in the reference label map L 1 (see Figure2)

  Again, we consider instead Ř = [b b b 1 | b b b 2 ], and prove that for any ε > 0 there exists e e e, e e e ∈ Ǧ, linearly independent, such that e e e < ε, e e e < ε. Let {.} stand for the fractional part function. We study the images of {Zb b b 1 } ∈ [-1/2, 1/2) 2 , where b b b 1 = a/c b/c denotes the first column of Ř. If b b b 1 contains irrational elements, the {Zb b b 1 } contains infinitely many distinct points. By compactness of [-1/2, 1/2] 2 , we can extract a subsequence ({n k b b b 1 }) k∈N , converging to some point in [-1/2, 1/2] 2 . Thus {(n k+1n k )b b b 1 } = {n k+1 b b b 1 } -{n k b b b 1 } converges to (0, 0). In particular, we can find integers m, p, q, where m = n k+1n k for k large enough, such that e e e = mb b b 1 + (p, q) ∈ Ǧ has norm smaller than ε 3 . Note now that the second column of Ř satisfies b b b 2 = 0 -1 1 -1 b b b 1 . Then we claim that e e e = 0 -1 1 -1 e e e = mb b b 2 + (-

Corollary 10 .

 10 If cos θ = 2a-b 2c and sin θ = √ 3b 2c where (a, b, c) is a primitive Eisenstein triple, the group Ḡ = G/Λ is cyclic and | Ḡ| = c.

  Fig.5. Comparison of the loss of information induced by digitized rigid motions on the hexagonal and the square grids. The red and blue curves correspond to the ratios between the areas of noninjective zones and the remainder range for the square and the hexagonal grids, respectively. The "x" (resp. "o") markers correspond to the values of the information loss rate 1-|U(S Z )| |S Z | (resp. 1 -|U(S L )| |S L | )

Fig. 6 .

 6 Fig.6. All the neighborhood motion maps κ∈Λ G U r (κ) for any angle in (α 0 , α 1 ) (top), (α 1 , α 2 ) (bottom, left), and (α 2 , α 3 ) (bottom, right), visualized by the label map L U 1 (see Figure2). Neighborhood motions maps are considered as graph vertices and linked by edges with respect to the adjacency relations of respective frames, i.e. each edge in the graph represents a side shared by two frames in the remainder range. The neighborhood motion maps which correspond to noninjective zones are surrounded by pink ellipses. The elements which are changed with respect to the rotation angles in the bottom figures are surrounded by black squares, while those which are not changed are faded. A high resolution version of the figure can be found via the URL: https://doi.org/10.5281/zenodo.820911

Similar attempt was made in[START_REF] Thibault | Rotations in 2D and 3D discrete spaces[END_REF], but with a different approach of measuring.