Jean-Charles Papin
email: jean-charles.papin@ens-cachan.fr

Christophe Denoual

• Laurent

Keywords: Simulation, dynamic load-balancing, graph partitioning, tasks, many-core, pair potential

Many applications of physics modeling use regular meshes on which computations of highly variable cost over time can occur. Distributing the underlying cells over manycore architectures is a critical load balancing step that should be performed the less frequently possible. Graph partitioning tools are known to be very effective for such problems, but they exhibit scalability problems as the number of cores and the number of cells increase. We introduce a dynamic task scheduling and mesh partitioning approach inspired by physical particle interactions. Our method virtually moves cores over a 2D/3D mesh of tasks and uses a Voronoi domain decomposition to balance workload. Displacements of cores are the result of force computations using a carefully chosen pair potential. We evaluate our method against graph partitioning tools and existing task schedulers with a representative physical application, and demonstrate the relevance of our approach.

Many physics simulation applications rely on large meshes where each cell contains a few elementary computing elements (e.g. particles, finite elements) and is linked to its neighboring cells. Figure 1 illustrates the concept of mesh and neighboring dependencies. The computing cost of each cell varies with the modeled material and the applied action (e.g., shock wave or distortion). Thus, distributing cells among computing units must both preserve locality of neighbors and dynamically balance computing load. Mesh partitioning can either be achieved by using generic task scheduling or graph partitioning tools.

Along with specific libraries [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF][START_REF] Kukanov | The foundations for scalable multi-core software in intel threading building blocks[END_REF]19] and languages extensions [START_REF]Openmp application program interface 3.0. Technical report[END_REF], many languages now integrate task support in their standard, such as C++11 with its future variables concept. Tasks are generally associated with a data set and can exchange data with other tasks, allowing to define affinity criteria for a given task. Therefore, a particular attention shall be paid to the place (i.e., the computing unit) where the task will be executed. The still increasing coreson-chip number [START_REF] Hemmert | Green hpc: From nice to necessity[END_REF] forces applications to spawn a large number of tasks so as to prevent cores from entering a starving state. This eventually leads to an important complexity as task scheduling is known to be an NP-Complete [START_REF] Ullman | Np-complete scheduling problems[END_REF] problem. Moreover, existing dynamic task scheduling policies minimize the overall computing time and/or data displacements. By doing so, they do not necessary ensure the compactness of the set of cells to be exchanged, thus possibly increasing the cost of communications when physical models involve many neighbor-to-neighbor communications.

On another hand, graph partitioning tools can be used to build compact sets of mesh. Several well known libraries exist [START_REF] Pellegrini | Scotch and libScotch 5.1 User's Guide[END_REF][START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF][START_REF] Devine | Design of dynamic load-balancing tools for parallel applications[END_REF] and provide highquality partitioning algorithms. The NP-Complete problem is handled by using multi-levels algorithms [START_REF] Hendrickson | A multi-level algorithm for partitioning graphs[END_REF] that successively reduce a graph into a smaller graph that can be partitioned with standard heuristics: the result is then successively propagated into the original graph. Such algorithms require to merge some nodes into bigger nodes, and thus, require to select graph nodes. This property is known to be limiting: in a distributed environment, such algorithms need to communicate a lot. With upcoming architectures (i.e., manycore) and with the still increasing computing clusters size, we definitely need algorithms that involve a limited amount of communications.

We present Spawn, a physical interaction inspired scheduler that produces compact and optimal Voronoi domains. In our case, Voronoi diagrams maximize per-core data locality by providing numerous advantages: cache usage improvements, more efficient NUMA-aware memory allocations/accesses and less Point-to-Point communications. This scheduler has the advantage of being efficient to compute and offers an automatic refinement. This scheduler can be used as a task scheduling algorithm as well as a mesh partitioning tool. This paper is organized as follows: section 3 presents the rationale behind using physical interactions for dynamic scheduling while section 4 focuses on implementation details and some preliminary evaluations. Following sections present a set of experiments that compare our approach with graph partitioning tools (section 6) and dynamic task schedulers (section 7). Concluding remarks and future work are discussed in the last section.

Context

Applications of physics modeling mainly encounter two load variation scenarios: randomly diffused load or a continuous load variation, both over time. We build a grid of tasks that corresponds to the modeled material with an analogy between sets of elements (atoms, finites elements) and computing tasks. In this case a task is essentially a dataset on which we regularly apply some computing functions. Figure 2 represents the set of cells on a multi-dimensional grid. Each cell has its computing cost of its own that can evolve over time. Distributing all these cells among computing units is complex and several tools or methods exist. As a consequence of neighboring dependencies (see fig. 1), we must distribute cells into compact sets. These sets must be as conservative as possible among successive scheduling steps in order to prevent data displacements. A common way to achieve this is to move sets boundaries (see fig. 3). However, such decomposition induces an unstable number

t = 0 t = 1
Fig. 3 Regular domain decomposition at t = 0 (left) and t = 1 (right). When moving upper and/or lower boundaries, the number of connected neighbors can strongly vary.

of neighbors over successive scheduling: moving upper and lower boundaries modify the number of connected neighbors too. Graph partitioning tools [START_REF] Pellegrini | Scotch and libScotch 5.1 User's Guide[END_REF][START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF] can also be used to distribute tasks. By exploiting the fact that tasks have strong affinities with their neighborhood (see fig. 1), such meshes can directly be represented by graphs where each node is a task, and where a link is a neighboring connection. Partitioning tools provide very efficient algorithms for very large regular and non-regular graphs, which can take into account affinity characteristics. Nevertheless, such tools do not perform well when dynamic scheduling is required. If graph refinement (i.e. the ability to update current partitions) is not supported, slight task weight variations usually lead to the generation of a completely different graph partition, which incurs numerous data transfers to re-assign tasks. Parallel algorithms involved into distributed versions of partitioning libraries [START_REF] Chevalier | PT-Scotch: A tool for efficient parallel graph ordering[END_REF][START_REF] Karypis | Parallel multilevel k-way partitioning scheme for irregular graphs[END_REF][START_REF] Devine | Design of dynamic load-balancing tools for parallel applications[END_REF], while offering refinement functionalities, suffer of scalability issues as we will see in section 6. In modern architectures (i.e. the Intel R Xeon Phi), the number of cores is significant (typically 61 cores for Intel R Xeon Phi 7100 series) and the number of available thread units doubles or quadruples this value. The amount of memory per core is consequently a limiting factor (less than 300 Megabytes that must be shared among threads) and low memory footprint load balancing mechanisms must be designed.

Load Balancing with Molecular Dynamics

One major contribution of this article is the integration of concepts coming from Particle Dynamics (PD) simulations, a type of N-Body simulations, into workflow scheduling. This section first describes how we use particles as Voronoi seed to tessellate the volume, and how an analogy with particles with electrostatics interactions is used to optimize the scheduling method. We will then explore how we target the optimal per-core load with this method and finally the way we translate the computing load into forces.

Using Particles to Assign Work to Processing Units

Set of Objects to Schedule

We define the elements to be scheduled between computing units as a grid of objects. Objects can be of different nature (atoms, finite element, computing task or graph node). The only two important properties of the grid are its dimensions (x, y, z, with x, y, z ∈ N) and its definition in the Euclidean space, allowing distance computations. Each cell of the grid has thus a specific coordinate (x, y, z, with x, y, z ∈ N) on the grid and a specific computing cost (see figure 4). Because of the geometric properties of the grid, neighboring cells of the grids can be seen as interconnected cells: a neighboring relation can be assimilated to a connection between two distinct tasks or graph nodes.

Xaxis Y -axis Fig. 4 The task grid model: the grid has a dimension (x, y, z, with x, y, z ∈ N) and is defined in the Euclidean space. Each cell has a specific coordinate and computing cost (right figure: darker cells have a more important computing cost than white cells).

We define a vCore as a virtual representation of a computing unit. When seen from the PD side, a vCore i is a particle: it has coordinates (in the Euclidean space, with x, y, z ∈ R), a velocity vector and force vector. A vCore can be positioned inside the previously-defined grid. When seen as a computing unit, a vCore 'i' has a computing load Q i defined as the sum of computing loads of its set of tasks. A relationship between electrostatics charges and the computing load is proposed in paragraph Particle dynamics (page 7) to complete the proposed analogy.

Tessellation

The assignment of cells to vCores (and thus, to computing units) resides in the partitioning of the previously defined grid. This work uses a Voronoï [START_REF] Okabe | Spatial tessellations : concepts and applications of Voronoï diagrams[END_REF] tessellation to assign cells to the underlying computing units. A Voronoï tessellation is a geometric partitioning based on distances between points. Voronoï diagrams gather a set of points around a central point: the site of a Voronoï cell (see figure 5). In our case, sites of Voronoï cells are defined by the location of our computing units (or vCore), and cells of the domain are our set of tasks (or node graph) to compute. Thus, distributing tasks to computing units is as easy as computing distances between set of points. The complexity of computing distances between points is linear in n (n • O(m), with n, the number of sites, and m, the number of cells), but specific methods can dramatically reduce this complexity to O(log(n)). Since the computing load Q i of a vCore i is the sum of tasks load inside its voronoï zones, obtaining load equilibrium between vCore therefore consists in finding the vCores locations (and thus the Voronï tessellation) that ensures equality between all Q i .

We chose to use Voronoï diagrams since they offer compact and convex sets and have geometric stability [START_REF] Reem | The geometric stability of voronoï diagrams with respect to small changes of the sites[END_REF]. They are used in many scientific applications (e.g., in physics [START_REF] Brostow | Construction of voronoï polyhedra[END_REF], geography [START_REF] Rhynsburger | Analytic delineation of thiessen polygons*[END_REF] and biology [START_REF] Blura | Biological shape and visual science[END_REF]). Interestingly, the compactness of the Voronoï sites arrangement plays an important role regarding the partitioned shapes, both in terms of surface of Voronoï zones and sites connectivity. For example, the Vornoï tessellation of a compact set of points in 2D leads to hexagonal shaped zones, and thus to the minimal perimeter/surface ratio compared to a simple square grid. In 3D, a compact set of points (e.g., points ordered following a face centered cubic lattice) tessellates the volume into regular rhombic dodecahedra, a space-filling polyhedra that also maximizes the surface to volume ratio.

Applied to our analogy between vCores and particles, a compact set of vCores minimizes the amount of surface, and therefore the number of cells located on the boundary of the vCore domain. Compared to a cubic arrangement of Voronoï sites, the number of neighbors of a vCore (i.e., the number of neighboring sites sharing a surface, a line or a point with a vCore) is reduced from 26 to 12, an interesting advantage when one want to limit the number of site-to-site communications.

The interaction between particles, presented in the following section, is designed to favor the compactness of the Voronoï sites arrangement as well as to allow for a handy control of the Q i optimizing procedure.

Particle Dynamics

A computing unit is associated to a Voronoï site (a particle) and to the corresponding tasks of the cells contained in its Voronoï zone. Particles arrangement thus define the way the objects to be scheduled are assigned. To design the interaction between particles, we have decided that the interactions should favor compactness of the site arrangement to minimize the number of site's connections, and that the calculation of the new partition should be done frequently so that its computational cost to achieve equilibrium must be marginal compared to the overall computations.

We consider in the following example of electrostatics interaction, a very interesting pair force candidate based on the Coulomb's interaction and producing compact structures (for homogeneous charges) at a very low computational cost. For two particles i and j, separated by a vector r ij = x i -x j , the electrostatic force is:

F ij = λ ij r ij |r ij | 3 , (1)
with λ ij a function of the loads (Q i and Q j ) of the interacting vCores to be defined in next section (note that λ ij is the product of particle's charge in a standard Coulomb's interaction).

Focusing on the equilibrium only, we use a damped relaxation of the forces by expressing velocities as a factor of applied forces:

v i = dx i dt = -αF ij (2)
with α, a positive scalar. We also lump α with pseudo-time increment into a simple scalar k:

x i (t + ∆t) = x i (t) -k j F ij ( 3
)
We rescale k for every step so that the distance x(t) -x(t -∆t) is a fraction of the cell box dimension, which ensures a rapid convergence to stable states.

Converging Toward Optimal Load Balance

As explained in previous section, we want to make vCores move so that they reach a position at which, once the tessellation step achieved, equilibrates the load between vCores. Besides the fact that this requires a specific force equation definition, one can note that the exact solution of this problem relies on a global knowledge: we compute forces between pairs of vCores, and we target an optimal load distribution that satisfies the load of all vCores. We define the optimal vCore load m as the sum of the load of all vCores, which is by definition the cost of the whole cell grid, and we divide it by N , the number of vCores. Therefore, m is the average computing load of the whole computing domain:

m = 1 N N i=0 Q i (4)
Since this optimal per-vCore definition would result in a poorly scalable centralized implementation due to need of knowing the load of the cells of the whole domain, we define an approximate per-vcore local optimal that only uses the load of surrounding cells:

m approx = 1 N s Ns i=0 Q i (5)
While this definition is close to the previous one (see equation 4), the resulting solution is no longer exact, and the convergence is toward a local optimal. Now that our targeted convergence state is defined, we need to define the force equation that takes into account our optimal load distribution.

Force Equation

Forces between particles are proportional to the interaction amplitude λ ij (see equation 1)), designed in this section so that the vCore load Q i tend to the average load m, targeted by all vCores. The difference m -Q i represents the distance between the load of a particle i and the average (optimal) load m. For two interacting VCores i and j, we take a force proportional to the sum of two distances (m -Q i and m -Q j ) instead of the product of charges C i C j in Coulomb's force definition:

F ij = λ r ij |r ij | 3 , with λ = 1 - Q i + Q j 2m . ( 6
)
This force definition, referred to as the MediumLoad potential, produces three kinds of forces: positive, negative and null forces. For two vCores i and j, when the Q i + Q j sum is smaller than the 2m term, the whole fraction Qi+Qj 2m is smaller than 1, and the whole λ is positive, leading to positive forces. Similarly, when the Q i + Q j sum is larger than the 2m term, the whole fraction is larger Fig. 6 Left: dipoles formation when using our force equation in its original definition (see equation 6. Right: interactions when using the improved equation version (see equation 7.) than 1, thus, the λ term is negative leading to negative forces. Lastly, when the Q i + Q j sum is equal to the 2m term, the λ is null, and forces are null. An equilibrium is thus achieved when all the loads Q i are equal to the target load m.

With two particles, the Q i + Q j sum is smaller than the 2m term only when the two vCores are under-loaded, and is larger than the 2m term only when the two vCores are over-loaded. When an under-loaded vCore repulses neighboring vCores through repulsive forces, its own Voronoï domain grows by collecting cells. In the same way, when an over-loaded vCore attracts its neighbors through attractive forces, other vCores grab cells from it.

Fixing Dipole Formation

The proposed force equation allows two particles to collapse into a dipole as long as their load are equal to the average load (see figure 6). This brings two major issues: vCores spinning and, when the distance tends towards zero, to a singularity (division by zero). To eliminate this degenerate state, we have then added a repulsive term, also known as a short repulsive term, ensures a minimal distance between two vCores:

F ij = λ r ij |r ij | 3 + ν r ij |r ij | 5 , with λ = 1 - Q i + Q j 2m . ( 7
)
with ν a parameter tuned to have a minimal distance between particles about the size of cell box dimension.

Implementation Details & Preliminary Evolution

Global Vs Local Targeted Load

As previously described, we have defined two ways to regulate the load between computing units. A global one, that requires the global knowledge of costs of all tasks of the domain and that requires to compute forces between all pair of vCores. The second one considers only local interactions, and local information. Thus, it relies on a neighbor list, and does the whole load balancing step by using only this local information (reduced set of cells, reduced set of vCores). Figure 7 presents these two force schemes. While this second option has good parallel opportunities (i.e., only local node-to-nodes communications), this scheme implies two drawbacks: we consider only the local particle system instead of N-Body system, leading to a non-perfect solution to initial load balancing problem. Therefore, the second disadvantage is that we could have a higher number of iterations of the particle dynamics main loop algorithm to reach the targeted optimal load.

Tasks Assignment: Centralized Vs Distributed

Task assignment is also a critical point because it depends on the position of vCores, i.e., we partition the cell domain by using the final position of vCores.

We use two distinct algorithms to perform cell assignment: a fully sequential one and a fully distributed one based on the only-node-to-node communication model. The first one (see algorithm 1) is to be computed by a single core, and thus, in a distributed memory environment it requires a lot of All-To-One communications. The second one, presented in details in section 4. is able to compute the Voronoï tessellation in parallel (i.e., each computing unit computes its own Voronoï tessellation) by using only communications to neighboring nodes. To do so, this algorithm needs to know the position of the local vCore inside the domain and the position of its neighbors. With these two information, a vCore can then evaluate whether a cell is inside its local Voronoï zone or not. We recall that this step consists in closest-element search. We don't want each vCore to iterate over all cells of the domain, and by using the bounding box of the Voronoï zone n, we reduce the set of cells on which we iterate to deduce the Voronoï zone n +1 .

Parallel Force Computation

Because of the weak impact of long-range interactions, we can avoid to compute forces between all pair of particles. Therefore, in parallel, we consider only first and second-rank neighboring interactions: every node knows the computing load of its neighboring nodes and can thus deduce all locally-applied forces by its neighbors. However, because of the of our potential definitions (see section 4.1), this parallel force computation step is limited by the way we deduce the optimal per-core load (i.e., fully centralized or with local interactions only).

Different Versions

In order to clearly evaluate our method, we have designed several versions of our load balancing model.

Version 2 (V2): distributed algorithms, exact solution

In this version, we use a distributed cell assignment and force computation, but with the global optimal load. The main convergence loop is presented in algorithm 3: In this version, the computeLocalF orce() step relies on a All-To-All communication in order to retrieve the global optimal m value presented in equation 6Version 3 (V3): distributed algorithms, local solution This version is exactly the same as the V2 but the computeLocalF orce() step no longer relies on a global communication scheme. Instead, we deduce the local optimal m of equation 6 by using the load information of our neighbors.

Algorithm 3 PSpawn Main Convergence Loop. We first need to synchronize (line 2) neighbor information (due to the newest local computing load). Local force and thus, new local position can be computed (lines 3 and 4). We then need to synchronize (line 5) position information before computing the local Voronoï domain (line 6). The last step (line 7) exchanges costs of tasks that leave/arrive into the local domain.

Benchmark Application

This section presents results of our scheduling library (Spawn). In order to evaluate our scheduler, we have developed a representative stencil-like model applications (see table 1). This benchmark application uses a real cell domain (with evolutive load, see figure 8) and provides an intensive memory usage with data exchanges between cells. Every cell has 4 direct neighbors (four cardinal directions), and requires data (a matrix of double) from its neighbors to compute its own matrix. These exchanged matrices are the result of the previous iteration. In a distributed memory (see section 6.1.1), this allows us to evaluate the number of MPI communications induced after the new domain decomposition and during computation (data exchange between meshes). We evaluate our scheduler with sequential partitioning tools and parallel partitioning tools (with partitioning refinement enabled). In section 7, we evaluate the same application, but in a shared memory, by comparing our scheduler against common task scheduling strategies. Our analysis focuses on cache miss rates. In a shared memory environment, memory accesses must be handled with care to avoid performance bottleneck. Due to hierarchical memory organisation (several NUMA nodes, several cache levels), the data locality, usage and reusage are mandatory to get good performance. In this section, we compare performance and data cache misses ratio and show how the Spawn library improves this ratio in both static and dynamic load evolution (see figure 8).

Spawn-V1 VS Graph Partitioning Tools

Figure 9 compares performance of Spawn (version V1) in relation to Scotch [START_REF] Pellegrini | Scotch and libScotch 5.1 User's Guide[END_REF] and Metis [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF] libraries. Scotch supports graph refinement and all experiments featuring Scotch exploit this functionality. Left figures refer to the static load case while right charts refer to the dynamic load evolution case, and, upper charts refer to the speedup while the bottom ones refer to the average number of L 1 data cache misses per thread for one iteration. All partitioning strategies produce compact sets, minimizing data transfers between tasks, but Spawn has better performance in both static and dynamic load variations. Graph partitioning tools may not support graph refinement (e.g. Metis) and thus, produce a new and totally different task distribution after a task cost variation. Due to hardware and Operating System behavior, a task cost may vary over time, even if it's the same task and it does the same type of computation. Hence, such partitioning tools produce different task distributions after each iteration. This involves data displacements between sockets and higher cache miss rates (see bottom charts in figure 9).

Distributed Memory: Application #2b

Distributed memory implies to send data over a transfer bus (a network, a Pci-E bus, and so forth). Therefore, high contentions can occur upon a large data transfer, and especially when several computing nodes exchange data at the same time. In this section, we evaluate Spawn V1 performance against Scotch is used with refinement strategies (SCOTCH graphRepart routine). Task scheduling or graph partitioning is achieved after each computing iteration. The Spawn cache miss ratio is used as a reference on bottom charts. sequential and distributed graph partitioning tools. We particularly focus on speedup and on the number of exchanged cells.

Sequential Partitioning Tools

Figures 10 and 11 present the same evaluation but in a distributed memory environment by using MPI libraries. Measurements focus on speedup and data transfer rates between each iteration. We use matrices of 128x128 doubles in a domain of 512x512 tasks. Tests run on a cluster of Intel R Xeon X7650 nodes, interconnected with an Infiniband QDR network. Figure 10 shows the average number of exchanged tasks on the MPI network (in percentage of the entire domain) and the resulting application speedup. Metis involves an important data exchange since it is not designed for dynamic partitioning: nearly the entire task domain is exchanged for each iteration while less than 5% of the domain is exchanged with Scotch and Spawn.

Distributed Partitioning Tools

Tools like PTScotch [START_REF] Chevalier | PT-Scotch: A tool for efficient parallel graph ordering[END_REF] and ParMETIS [START_REF] Karypis | Parallel multilevel k-way partitioning scheme for irregular graphs, department of computer science[END_REF] account the actual partition to produce a new partition with the newest task costs. Scotch does not. Zoltan [START_REF] Devine | Design of dynamic load-balancing tools for parallel applications[END_REF] is another partitioning tool that can use not only its internal algorithms but also Scotch/PTScotch or ParMETIS to compute a graph partitioning. One limitation of these tools is that they need the identifier of each node that owns every local tasks and every neighboring local or remote tasks. This implies to maintain on each MPI node, and after every calls to a partitioning routine, the list of neighboring nodes, and the list of tasks they share.

The next evaluation focuses on Zoltan and ParMETIS (with Zoltan) tools. We evaluate Zoltan with PHG and RCB partitioning algorithms. PHG is the internal Parallel Hypergraph and Graph partitioning method. It supports initial partitioning, refinement and re-partitioning operations. Both repartitioning and refinement operations reuse the current partition, but refinement operation is stricter. With the RCB method we use one internal algorithm option, RCB REUSE, that indicates whether previous cuts should be used as initial guess for the current RCB partition. Figure 11 shows data transfer rates and speedups of Zoltan (with PHG, RCB and ParMETIS) in comparison with Spawn. We can see that Spawn performance is similar to graph partitioning and geometric algorithms when the network is not used (i.e. one computing node used). For a few number of MPI nodes, Spawn is still better than Zoltan or ParMETIS, while, for a higher number of MPI nodes, RCB algorithm outperforms other algorithms. More precisely, if we compare the RCB algorithm to Spawn, we can see (left chart), that Spawn produces less data transfers (1-3% of the entire domain versus 5-35% for RCB). This explains that our algorithm involves a better speedup when the network is not used (time spent on internal communication is regained by the time not spent on exchanging data). On the other hand, with a higher number of nodes, the centralized behavior of Spawn, i.e., the knowledge of the average load and its single-thread implementation, is penalizing and collective MPI operations considerably slow down the global computation. This observation is however incomplete if we compare Spawn to the PHG method or to ParMETIS. In comparison with Spawn, data transfer rates are more important for PHG, and really similar to ParMETIS. Thus, we could expect better results for, at least, ParMETIS. Considering that partitioning tools have scalability limitations because of internal synchronizations, exacerbated by the network, we can consider that data transfers are less penalizing than scalability limitations, which explains PHG and ParMETIS results in figure 11.

Evaluation of the Distributed Algorithm

This section presents and analyzes induced performance in terms of speedup and number of communications of the Spawn library by using version 2 (V2) and 3 (V3). Experiments are achieved on a cluster of 32 Intel R Xeon Xeon E5-2698v3 nodes (Haswell architecture, 2x16 cores, 128 Gio of RAM, for a total of 1,024 cores), interconnected with an Infiniband QDR network. We use the same application set as presented in table 1.

We compare performance of the Spawn library in relation with other parallel graph partitioning tools. As in previous section we compare performance against the Zoltan library. Two internal and one external algorithms are used: RCB (the geometric graph partitioning algorithm), PHG (the graph partitioning algorithm), and the external ParMetis library. A test configuration consists in a grid of 512x512 cells. Each cell contains a matrix of 128x128 double values (128Kio per matrix). Figures 12 and13 show the cost of our iterative model. The first one presents speedup of our benchmark application over an increasing number of MPI nodes and with an increasing number of PSpawn convergence iterations (I parameter). The second one presents the time spent in main parts of our distributed algorithm. We can first see that the application speedup, while greatly enhanced with our PSpawn algorithm, strongly depends on the number of iterations: we lose nearly 35% speedup points with 40 iterations regarding the one with one single iteration. We can also notice that we have a brutal performance drop when using 512 (and more) MPI nodes. One can expect better performance when increasing the number of iterations of the model (I parameter) but results show that best results are achieved with I = 1. This is due to the way our benchmark application is designed: it has a few computing ratio regarding its communication volume (each cell requires data from its 4 direct neighbors) and it computes a new cells distribution between each of its computing iteration. Therefore, the scheduling step is limiting. In order to confirm this point, figure 14 presents a speedup over 1,024 MPI nodes with the same domain size (512x512 cells) but with a more important amount of computation (x5 factor). We can clearly see that the scheduling step is less penalizing than before. Fig. [START_REF]Openmp application program interface 3.0. Technical report[END_REF] PSpawn performance with a more CPU-bound application (x5 regarding other tests). In this case, scheduling is less penalizing.

Figure 13 shows that most part of the computing cost (more than 80%) of our algorithm remains in the Voronoï tessellation, and since we need to compute a Voronoï tessellation for every single iteration, this explains why the speedup is penalized when we increase the number of iterations. We also see that, as we increase the number of MPI nodes, the part of communications is more important (for both force neighbor communications and force computations). On Figure 15, which shows real time values, we see (left chart) that the Voronoï tessellation cost dramatically decreases with the increase of the number of MPI nodes. The time spent in neighbor-to-neighbor communication (middle chart) slightly increases for a small number of processes (between 1 and 16) but remains however constant for a higher number of MPI nodes. This is a consequence of our peer-to-peer model: when using a higher number of nodes, communications are limited by the nearly constant number of neighbors. We can also compare the impact of the centralized property of the V1 potential: it requires, during the force computation step, to know the whole task cost of the cell domain, leading to an All-to-all MPI communication. Even if figures 13 and 15 show that we clearly reduce the communication cost during the force computation with the V2 potential, figure 12 shows that removing the global communication required by the V1 potential does not have a significant impact on overall performance. Deeper experiments have shown that the V2 potential is unstable (non-constant quality) and induces several cell movements between each application computing iteration (see figure 16).

Finally, figure 16 presents the average number of exchanged tasks after the computation of new partition. As one can notice, this number slightly increases within the number of iterations. This is an expected behavior since vCores displacement is related to the number of iterations. This number is however really similar to the one induced by using graph and geometric partitioning tools. We can also notice the instability of the V2 potential: since vCores do not target a global computing load, they have independent local displacements, leading to oscillations during displacements. In the end, these oscillations induce a lot of non-desirable cells, and therefore data displacements between computing nodes.

Comparison with Affinity-Based Scheduling Policies

We have extended the StarPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] runtime with our scheduler to evaluate the quality of the Spawn (V1) performance against common task scheduling policies. Data block of tasks are managed by StarPU through the specific starpu block data register routine, which enables automatic and backgrounded data transfers. We evaluate our scheduler on a two-socket Intel R Xeon E5-2680 (IvyBridge, 2x10 cores, 2.8GHz) machine. We did a hundred of runs with a domain of 100x100 tasks, where each task contains a matrix of 8x8 doubles. We evaluate performance in both static and dynamic load variations over time.

Integration

StarPU comes with a variety of scheduling policies, but also offers the capability to define new ones. This permits to manage StarPU workers (abstract representation of computing units) directly. Spawn is used by StarPU through an intermediate meta-scheduler that implements interfaces required by StarPU (the starpu sched policy C-structure). When tasks are submitted to StarPU, they are forwarded to the scheduling layer, and thus, to our meta-scheduler. This meta-scheduler asks our library to which worker this task should be assigned (more precisely, to which Voronoï domain this task is attached). Regular grids or random positions are used to initialize the vCore set, and therefore the first task distribution. Once tasks are distributed among StarPU workers, their execution (e.g. the computing time) is monitored and then sent to our library. The scheduling is computed with nearly no additional cost (see figure 17): it is computed by the first inactive StarPU worker, and as long as another worker is still working. Thus, we only compute a task distribution when an imbalance exists. This has also the advantage of allowing us to set P = 1 and I = 1 as parameters of our algorithm (which can be interpreted as "as long as an imbalance exists, do one iteration to reach the most optimal distribution"). As explained, a new task distribution is based on the task information (i.e. computing cost) of the previous task execution. Therefore, the first two load balancing steps are not optimal.

Results

Figure 18 compares performance of Spawn (V1) in relation to predefined StarPU schedulers. Left charts refer to the static load case while right charts refer to the dynamic load case. Upper charts refer to the speedup while the bottom ones refer to the ratio of L 1 data cache misses per thread (per iteration) in comparison with Spawn. We use three internal StarPU schedulers: Eager, DM and DMDA since other schedulers are only useful with heterogeneous machines. Eager uses a simple FIFO-based greedy policy to schedule tasks. DM and DMDA use performance models allowing them to perform smart task placement, with the objective to minimize the overall execution time. DMDA works like DM, but takes into account data transfer time during task placement. More information about DMDA can be found in the StarPU HandBook [START_REF]INRIA Runtime. Starpu handbook[END_REF].

In both static and dynamic load variation cases, Spawn is close to StarPU schedulers performance and outperforms them with the growth of the number of threads. Bottom parts of figures show that Spawn dramatically reduces the number of L 1 cache misses with a factor between two and four. This is due to the Spawn ability to provide geometrically compact sets of tasks, reducing data transfers between threads. On the other side, if we compare cache misses with Eager, DM and DMDA, we note that Eager doubles the number of cache misses that DM and DMDA induce, but has better performance. This is related to the behavior of these schedulers: by using a task queue, Eager naturally maximizes StarPU workers activity, in opposition to DM and DMDA that distribute tasks in advance as soon as they are ready in order to minimize overall execution time. Since our application uses intensive memory accesses, the execution time is strongly related to memory accesses, and task execution time depends on the location where tasks are executed. Therefore, assigned tasks by DM or DMDA constitute a set of non-related tasks, like they would have been assigned with Eager. Hence, with DM or DMDA, StarPU workers are imbalanced, since the execution time of assigned tasks does not reflect their real execution time. A deeper analysis of StarPU execution trace (FxT [START_REF]Fxt library for execution traces generation[END_REF] traces, see fig 19) confirms this point. With this information, it appears that common task schedulers are unable to perform efficiently with such task configurations (neighboring and data dependencies).

Conclusions

We introduce a new approach to dynamically schedule tasks within domaindecomposition based applications, that borrows many ideas from molecular dynamics simulations. We combine physical particle interactions and Voronoi domain decomposition to virtually move cores over a 2D/3D mesh to partition the set of application tasks. We have designed three versions with either centralized or distributed properties. In a shared memory environment, or with a few number of MPI nodes (around 256), experiments show that the fully centralized version outperforms well known graph partitioning tools and common task schedulers. Versions with distributed algorithms improve scalability while keeping important properties (data locality, automatic and cheap partition refinement) but slightly reduce the scheduling precision. Experiments show that our approach preserves tasks locality and remains stable in the presence of significant load variations. Our approach compares favorably to parallel graph partitioning tools (Zoltan, ParMetis) by introducing less data transfers and a faster refinement step. In the near future, we plan to introduce new potentials that would handle heterogeneous nodes with either hybrid Cpus (i.e. Arm big.Little architecture), with mixed computing units (Gpus + Cpus), or with different Cpus architectures. We also aim at supporting more general input meshes, in order to extend the scope of our library to applications relying on irregular meshes.

Fig. 2

 2 Fig.2Grid of tasks with an evolving task cost. Each square represents a task, and the brightness represents the cost of the task: dark cells are more costly than white ones. Each task contains a set of the physics elements of the simulated domain.

Fig. 5

 5 Fig. 5 Voronoï tessellation. On the left figure, Voronoï sites are positioned on the cell domain. The right figure shows the cell distribution into Voronoï cells, characterizing the task distribution over underlying computing units.

Fig. 7

 7 Fig.7Two interaction schemes. Left: fully centralized one based on global knowledge. In this scheme, every vCore needs to know the cost of all cells, and to communicates with all others vCores. Right: local interactions only. Here we consider only neighboring nodes and cells.

Version 1 (Algorithm 2

 12 V1): sequential algorithms, exact solution In this version, we use only a single computing unit to compute the force iterations and the cell distributions between all other computing units. As presented in algorithm 2, we have to perform global All-To-One communication steps (blocking step). Voronoï Tessellation ... receiveCellCostInformationFromALL () mainConvergenceLoop () // load balancing step for v in setOfvCores do sendCellAssignment () end for ...

Fig. 8

 8 Fig.[START_REF] Hendrickson | A multi-level algorithm for partitioning graphs[END_REF] Generated load evolution: each square represents a task, and the color represents the cost of the task: dark cells are more costly than the white ones, i.e, we do more computations on internal matrices.

Fig. 9

 9 Fig.[START_REF] Karypis | Parallel multilevel k-way partitioning scheme for irregular graphs, department of computer science[END_REF] Resulting speedups (top) and L 1 data cache misses (bottom) with graph partitioning tools (Scotch, Metis) and with Spawn on an Intel R IvyBridge Xeon E5-2680 computing node. Scotch is used with refinement strategies (SCOTCH graphRepart routine). Task scheduling or graph partitioning is achieved after each computing iteration. The Spawn cache miss ratio is used as a reference on bottom charts.

Fig. 11

 11 Fig. 11 Resulting data transfer rates (left) and speedups (right) with Zoltan (RCB, ParMetis (through Zoltan) and with Spawn. One can notice scalability issues of graph partitioning tools (PHG and ParMETIS).

Fig. 12

 12 Fig. 12 Cost of the PSpawn iterative model, successively increasing the number of iterations used to compute a partition. Top left: 1 iteration. Top right: 10 iterations. Bottom left: 20 iterations. Bottom right. 40 iterations.

Fig. 13

 13 Fig.[START_REF] Okabe | Spatial tessellations : concepts and applications of Voronoï diagrams[END_REF] Cost repartition in our distributed algorithm by using our two distinct version for different number of MPI processes: the V1 (top) and the V2 (bottom). Neighboring communications consist in exchanging local node information (load, position) and in the cell cost exchange between each local Voronoï tessellation: we need to send (or receive) cost of cells that leave (or arrive in) the local domain.

Fig. 15

 15 Fig.[START_REF] Pellegrini | Scotch and libScotch 5.1 User's Guide[END_REF] Cost evolution of the Voronoï tessellation, the communication ratio and the force computation steps of the distributed algorithm. Note the 10 5 scale factor for the cost of the Voronoï Tessellation chart (right). The vertical red line shows the limit after which the network is involved. Dashed lines represent the cost evolution when using the V2 algorithm.

Fig. 16

 16 Fig. 16 Comparison of average induced communication per MPI node after the computation of a new cell distribution with a variable number of iteration steps (I parameter). Top left: 1 iteration. Top right: 10 iterations. Bottom left: 20 iterations. Bottom right. 40 iterations.

Time

Fig. 17

 17 Fig.[START_REF] Rhynsburger | Analytic delineation of thiessen polygons*[END_REF] Timeline representing StarPU calls (grey), task computation (green) and partitioning computation with Spawn (red). The scheduling time is fully hidden by computing imbalance.

  Fig.[START_REF] Rhynsburger | Analytic delineation of thiessen polygons*[END_REF] Timeline representing StarPU calls (grey), task computation (green) and partitioning computation with Spawn (red). The scheduling time is fully hidden by computing imbalance.

Fig. 18

 18 Fig.[START_REF]INRIA Runtime. Starpu handbook[END_REF] Resulting speedups (top) and L data cache misses (bottom) with internal StarPU schedulers (Eager, DM, DMDA) and with Spawn on an IvyBridge Xeon E5-2680 computing node. The performance drop for nearly 10 core indicates the NUMA-effect. With Spawn, this effect is delayed regarding other scheduling strategies. The Spawn cache miss ratio is used as base on bottom charts.

Fig. 19

 19 Fig. 19 FxT traces with the DM scheduler. Green elements refer to task execution while red elements refer to inactivity of StarPU workers. One row symbolises the activity of one StarPU worker.

Table 1

 1 List of benchmark applications and their properties

	App No Parallelism Matrix size	Memory	Focus
	#a	pThreads	8x8	Shared	Cache misses
	#b	MPI	128x128	Distributed Data transfers

  can perform parallel graph partitioning. ParMETIS supports refinement graph partitioning, i.e., it takes into Metis) and with Spawn. Scotch is used with refinement capabilities. Scheduling or partitioning is achieved after each iteration.

						Spawn	Scotch	Metis
		100
	Exchanged tasks (%)	10				Speedup	10 30 20
		1
							0
		0	200	400	600	800 1,000	0	200	400	600	800 1,000
			Number of MPI processes		Number of MPI processes
	Fig. 10 Resulting data transfer rates (left) and speedups (right) with graph partitioning
	tools (Scotch,


Simone De Marco

Minh-Duc Hua

Robert Mahony

Tarek Hamel

Homography Estimation of a Moving Planar Scene from Direct Point Correspondence

Keywords: Homography, Internal Model Observer, Nonlinear Observer

Homography estimation is an important task in robotics applications such as landing on moving platforms, docking and refueling, building inspection, etc. Non-linear observers depend on an estimate of the group algebra velocity as a feedforward term to minimize lag in the filter response. When both the camera and the scene are moving, and for perspective image constructs such as the homography, it is often impossible to directly measure the required group velocity. This paper proposes a solution for the case where the motion is periodic, or approximately periodic, with known period. The approach is based on the internal model principle, where the internal model can be expanded to include a sufficient harmonics of the desired period in order to model complex periodic motions. The novelty of the work lies in formalizing the internal model for observer design to the non-compact Lie group SL(3) and providing a demonstration of the effectiveness of the approach with realworld examples.

I. INTRODUCTION

Recent years have witnessed an increase in the use of autonomous robotic vehicles for a wide range of civilian and commercial applications. There is a growing interest in applications where the autonomous vehicle physically interacts with the surrounding environment. These operations include landing on moving platforms, docking and refueling, building inspections and search and rescue missions. Obtaining a good measure of the position of the vehicle during such missions is a crucial to the system performance. Global Positioning Systems (GPS) are widely used for navigation and control in open areas, however, reliable and accurate GPS localization in the applications of interest is not possible due to tall buildings, urban canyons, dense material such as concrete and steel and many other factors that can block or heavily degrade the GPS signals. Moreover, many missions require localization with respect to some particular features in the environment instead of an absolute position. For these reasons, the study of alternative localization systems has become popular in the robotics and control community. Examples of such localization systems are 2D laser scanner (LIDAR), radar, cameras (monocular, stereo, RGB-D, event), etc. Active sensor systems such as LIDAR are well suited for large-scale vehicles, however since they are usually heavy, expensive and power hungry they S. de Marco and M.-D. Hua are with I3S, Université de Nice, Sophia-Antipolis, CNRS, France sdemarco(hua)@i3s.unice.fr R. Mahony is with the Australian Centre for Robotic Vision, Australian National University, Robert.Mahony@anu.edu.au.

T. Hamel is with I3S, Université Côte d'Azur, Institut Universitaire de France, CNRS, Sophia Antipolis, France, thamel@i3s.unice.fr. cannot be used on small-scale vehicles. In contrast, cameras are cheap, lightweight, small and convenient to mount and provide an information rich measurement. Therefore, visionbased systems remain one of the most promising technologies especially for small-scale unmanned aerial vehicles (UAV). Vision-based control algorithms are known in literature as visual servoing [START_REF] Chaumette | Visual servo control part I: Basic approaches[END_REF], [START_REF] Chaumette | Visual servo control, part II: Advanced approaches[END_REF], [START_REF] Malis | 2-1/2D visual servoing[END_REF].

Any two images of the same planar scene taken from two different camera poses are related by a projective mapping known as a homography. Homographies play a major role in rigid-body pose estimation [START_REF] Chen | Homographybased visual servo tracking control of a wheeled mobile robot[END_REF], [START_REF] Scaramuzza | Appearance-guided monocular omnidirectional visual odometry for outdoor ground vehicles[END_REF], unmanned vehicle navigation [START_REF] Nguyen | Inertial-aided homography-based visual servo control of autonomous underwater vehicles without linear velocity measurements[END_REF], [START_REF] Plinval | Visual servoing for underactuated VTOL UAVs: a linear, homography-based framework[END_REF], UAV autonomous approach and landing [START_REF] Gonc ¸alves | Homography-based visual servoing of an aircraft for automatic approach and landing[END_REF], [START_REF] Shakernia | Multiple view motion estimation and control for landing an unmanned aerial vehicle[END_REF], to name only a few applications. Using homographies is effective when the surrounding environment is composed mainly of planar surfaces for example in man-made environments, urban scenes and also when the relief of the scene viewed is negligible compared to the distance of the camera from the scene. Estimating homography from pairs of images has been extensively studied in the computer vision literature [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. Classical algorithms estimate homographies by minimizing a suitable cost function given a set of image feature (points, lines, conics) correspondences between image pairs. For such algorithms, there is no difference whether the images considered are generated by a temporal correlated sequence (such as in a video) or are taken independently. When a robotic vehicle moves, however, the camera captures a temporal sequence of images for which the homography varies continuously. This has lead to the development of a number of non-linear continuous-time observers [START_REF] Hamel | Homography estimation on the special linear group based on direct point correspondence[END_REF], [START_REF] Hua | Explicit complementary observer design on special linear group SL(3) for homography estimation using conic correspondences[END_REF], [START_REF] Hua | Point and line feature-based observer design on SL(3) for Homography estimation and its application to image stabilization[END_REF], [START_REF] Malis | Dynamic estimation of homography transformations on the special linear group for visual servo control[END_REF] that exploit the structure of the Special Linear group SL(3), a Lie-group isomorphic to the group of homographies [START_REF] Benhimane | Homography-based 2D visual tracking and servoing[END_REF]. A disadvantage of non-linear observers is that they rely on group velocity measurements to propagate the state estimate with time. Although robotic systems are usually provided with an inertial measurement unit, and other velocity measurements such as doppler GPS and optical flow, the homography velocity depends also on scene parameters and is difficult to reconstruct. Furthermore, if the scene itself is moving, such as for landing on a moving platform, ego motion estimation of the robot by itself cannot estimate the velocity of the relative homography. The only algorithms the authors are aware of that address this issue use an integral adaptive term [START_REF] Hamel | Homography estimation on the special linear group based on direct point correspondence[END_REF], [START_REF] Hua | Point and line feature-based observer design on SL(3) for Homography estimation and its application to image stabilization[END_REF], [START_REF] Malis | Dynamic estimation of homography transformations on the special linear group for visual servo control[END_REF] to estimate 'constant' unknown velocity. This assumption is restrictive since constant group velocity motions for the homography group are not intuitive motions for a rigid-body robotic system with a mounted camera moving over a moving scene, and it may introduce significant lag in the filter response when used in real-world scenarios [START_REF] Hamel | Homography estimation on the special linear group based on direct point correspondence[END_REF], [START_REF] Hua | Point and line feature-based observer design on SL(3) for Homography estimation and its application to image stabilization[END_REF], [START_REF] Malis | Dynamic estimation of homography transformations on the special linear group for visual servo control[END_REF].

The problem of regulating the output of a system in order to track an unknown reference trajectory while rejecting disturbances is known as an output regulation problem (or servomechanism problem) [START_REF] Davison | The robust control of a servomechanism problem for linear time-invariant multivariable systems[END_REF], [START_REF] Francis | The linear multivariable regulator problem[END_REF]. The case where the desired reference can be modeled as the solution to an ordinary differential equation with unknown initial condition can be solved by internal model-based regulators. This problem has been extensively studied over the last thirty years in both linear [START_REF] Davison | The robust control of a servomechanism problem for linear time-invariant multivariable systems[END_REF], [START_REF] Francis | The linear multivariable regulator problem[END_REF], [START_REF] Francis | The internal model principle of control theory[END_REF] and non-linear context [START_REF] Isidori | Output regulation for nonlinear systems[END_REF], [START_REF] Byrnes | Output regulation of Uncertain Nonlinear Systems[END_REF], [START_REF] Serrani | Semi-global nonlinear output regulation with adaptive internal model[END_REF], leading to the so-called Internal Model Principle. The central part of the regulator design incorporates an internal model unit, which is a dynamical autonomous system that generates in steady state the feed-forward input that keeps the tracking error identically to zero. Recently, the theory of non-linear observers has been shown to play a central role in the output regulation problem. In particular, adaptive observers [START_REF] Bastin | Stable adaptive observers for nonlinear time-varying systems[END_REF] and high-gain observers [START_REF] Gauthier | Deterministic Observation Theory and Applications[END_REF], [START_REF] Astolfi | A high-gain nonlinear observer with limited gain power[END_REF] have been successfully used for the design of adaptive internal models [START_REF] Delli Priscoli | A new approach to adaptive nonlinear regulation[END_REF] and robust internal models [START_REF] Isidori | Robust design of nonlinear internal models without adaptation[END_REF], [START_REF] Bin | Robust internal model design by nonlinear regression via low-power high-gain observers[END_REF].

In this paper, we take advantages of recent results [START_REF] Marco | Output regulation for systems on matrix lie-group[END_REF] on output regulation for systems on matrix Lie-groups to develop an internal model-based observer design for the homography estimation problem. The unknown group velocity, that acts like an unknown reference trajectory, is modeled as belonging to a family of trajectories generated as solutions of an autonomous linear system on the group Lie algebra. The proposed observer, in the same spirit of linear internal model principle, embeds a copy of the exogenous dynamics (posed on the Lie algebra) and incorporates an innovation to synchronize these dynamics with the exogenous plant providing a driving term for the velocity reference of a standard kinematic observer on the Lie group. The novelty of the paper, which is based on the direct point formulation proposed in [START_REF] Hamel | Homography estimation on the special linear group based on direct point correspondence[END_REF], is its ability to learn complex feed-forward velocity profiles to reduce lag in the homography tracking where motion is periodic. This has application in a wide range of applications such as landing an aerial robot on an oscillating platform or stabilizing periodic camera shake. A preliminary version of the theoretical results of this paper was presented in [START_REF] Marco | Homography estimation of a moving planar scene from direct point correspondence[END_REF]. Here we provide a more extensive stability analysis extending the local asymptotic results in [START_REF] Marco | Homography estimation of a moving planar scene from direct point correspondence[END_REF] to local uniform exponential results as well as providing a more detailed experimental study and full details of the derivations. The experimental results provided demonstrate excellent performance and robustness in presence of challenging conditions and outperform the state-of-the-art results documented in [START_REF] Hamel | Homography estimation on the special linear group based on direct point correspondence[END_REF].

The paper is organized into five sections including the present introduction and an appendix. Section II introduces the notation and formulation. The non-linear observer is detailed in Section III. Simulation results are provided in Section IV. Section V provides a detailed description of the software implementation of the observer. In the same section an application of the proposed approach for the tracking of an oscillating platform, with comparisons of the observer proposed in [START_REF] Hamel | Homography estimation on the special linear group based on direct point correspondence[END_REF], is presented. A video link of the experiment (https: //goo.gl/9Y3KsB) is provided as supplementary material. The Appendix contains proofs of some useful properties for the observer design, a rigorous analysis of the local exponential stability of the proposed observer, and an important lemma that is not in the direct line of argument for the main result.

II. PROBLEM FORMULATION

A. Notation and mathematical identities

The set of all m×n matrices whose entries are real numbers is denoted R m×n . For any A ∈ R n×n , det(A) and tr(A) denote the determinant and the trace of A respectively.

The Special Orthogonal group is denoted SO(3). The Lie algebra associated to the Special Orthogonal group, denoted by so(3), is the set of 3 × 3 skew symmetric matrices

so(3) := Ω × ∈ R 3×3 |Ω × + Ω × = 0 .
The Lie algebra so(3) with the matrix commutator (Lie bracket)

[•, •] is isomorphic to R 3 with the cross product. Let Ω = [Ω 1 , Ω 2 , Ω 3 ] ∈ R 3 , then the matrix Ω × = 0 -Ω3 Ω2 Ω3 0 -Ω1 -Ω2 Ω1 0
is the skew symmetric matrix associated to the cross product

Ω × v = Ω × v, for any v ∈ R 3 .
The Special Linear group SL(3) and its associated Lie algebra sl(3) are defined by

SL(3) := H ∈ R 3×3 | det(H) = 1 , sl(3) := U ∈ R 3×3 | tr(U ) = 0 .
For H ∈ SL(3) and U ∈ sl(3), the adjoint operator is a mapping Ad : SL(3) × sl(3) → sl(3) defined by

Ad H U := HU H -1 .
For any A, B ∈ R n×n , A, B = tr(A B) defines an inner product on R n×n , and the Frobenius norm • F is defined by

A F := A, A .
Let P sl3 (•) denote the unique orthogonal projection of R 3×3 onto sl(3) with respect to the trace inner product, one has

P sl3 (A) = A - 1 3 tr(A)I ∈ sl(3), ∀A ∈ R 3×3 ,
with I the identity matrix. Let ∧ denotes the mapping ∧ : R 8 → sl(3) that maps the vector v ∈ R 8 to an element of sl(3)

v ∧ := 8 j=1 v j B j where {B 1 , . . . , B 8 } is a basis of sl(3). For instance, denoting by b 1 = [1, 0, 0] , b 2 = [0, 1, 0] , b 3 = [0, 0, 1]
, the following eight generators form a basis for the lie algebra sl(3)

B 1 := b 1 b 1 -I 3 , B 4 = b 1 b 2 , B 7 = b 1 b 3 , B 2 := b 2 b 1 , B 5 = b 2 b 2 -I 3 , B 8 = b 2 b 3 , B 3 := b 3 b 1 , B 6 = b 3 b 2 . ( 1
)
The operator vec ∨ : sl(3) → R 8 denotes the inverse of the (•) ∧ operator, namely

vec ∨ (v ∧ ) = v, ∀v ∈ R 8 .
For any A ∈ R n×n , vect(A) ∈ R n 2 denotes the column vector obtained by the concatenation of columns of the matrix A as follows

vect(A) = [a 1,1 , • • • , a n,1 , a 1,2 , • • • , a n,2 , • • • , a 1,n , • • • , a n,n ] .
The matrix representation of the composition of the linear maps (vect •∧) : R 8 → sl(3) → R 9 is denoted by the full columns rank matrix J∧K ∈ R 9×8 . Let Q sl3 ∈ R 8×8 denote the symmetric positive matrix Q sl3 := J∧K J∧K.

For any p ∈ S 2 , the projection π p := (I -pp ) onto the tangent space of the unit sphere S 2 at point p satisfies π p p = 0.

B. Perspective Projection

Let B denote projective coordinates for the image plane of a camera, and {B} its frame of reference. Let ξ ∈ R 3 denote the position of the body-fixed frame {B} with respect to the inertial frame {I} expressed in {I}. The orientation of the frame {B} with respect to the inertial frame {I} is given by the rotation matrix R ∈ SO(3). Let T denote the projective coordinates of reference image of the target plane, and let {T } denote its right-hand frame of reference. Let (ξ T , R T ) ∈ R 3 × SO(3) denote the configuration of the frame {T } with respect to the inertial frame {I}.

The coordinates of a single point in the inertial frame P I ∈ {I}, in the camera frame P B ∈ {B} and in the target frame P ∈ {T } (see Figure 1) are related by

P I = R T P + ξ T , (2a)
P I = RP B + ξ. ( 2b
)
The orientation of the body-fixed frame {B} with respect to the target frame {T }, is denoted by the rotation matrix R : {B} → {T }, while the relative position of the frame {B} with respect to the target frame {T } expressed in the frame {T } is denoted by ξ, one verifies that

R = R T R, (3a) ξ = -R T (ξ T -ξ). (3b)
Substituting (2a) into (2b) one has

P B = R P -R ξ (4)
as relation between the coordinates of the same point in the target frame and in the body-fixed frame.

Considering the classical pinhole camera model [START_REF] Faugeras | Three-dimensional Computer Vision: A Geometric Viewpoint[END_REF], [START_REF] Ma | An Invitation to 3-D Vision: From Images to Geometric Models[END_REF], 3D points can be projected on the image planes B and T . Let p ∈ B denote the image of a point when the camera is aligned with the body-fixed frame {B} and let p ∈ T denote the image of the same point when the camera is aligned with the target frame {T }, then 1 p ∼ = K P , p ∼ = KP B 1 The notation " ∼ =" is defined in the homogeneous sense, i.e. up to a nonzero scalar factor [START_REF] Ma | An Invitation to 3-D Vision: From Images to Geometric Models[END_REF]. where K is the upper triangular camera matrix which depends on intrinsic parameters of the camera, namely focal length, pixels aspect ratio and principal point offset.

C. Homography matrix and the Special Linear Group SL(3)

Let η denote the unit normal vector pointing towards the target plane expressed in {T }, and let d denote the orthogonal distance of the plane to the origin of {T }.

Due to the fact that all target points Pi for i = {1, . . . , n} lie on a single planar surface one has

η Pi d = 1,
it follows from (4) that

P Bi = R -R ξ η d Pi , (5)
and thus, using the perspective projection, one obtains

pi ∼ = K R -R ξ η d -1 K -1 p i . (6)
The projective mapping G : B → T that maps pixel coordinates from B to T given by

G = ζK R -R ξ η d -1 K -1 ,
where ζ represents a scaling factor, is known in literature as homography matrix. Let η (resp. η I ) denote the normal to the target plane expressed in the frame {B} (resp. {I}), and let d (resp. d I ) denote the orthogonal distance of the planar surface to the origin of the body-fixed frame {B} (resp. inertial frame {I}). One verifies that:

η = R η I = R η, (7)
d = d I -η I ξ = d -η ξ. (8)
It follows that the homography matrix can be written as If the camera is calibrated the homography matrix can be scaled

G = ζK R + ξη d K -1 .
H = K -1 GK = ζ R + ξη d . (9)
This mapping is referred to as Euclidean homography since it maps Euclidean coordinates of the scene's points from {B} to {T }. Indeed, the Euclidean homography contains the relative pose information of the camera body-fixed frame {B} with respect to the target frame {T }. From now on we restrict our attention to the Euclidean homography H.

Since the homography H is only defined up to a scale factor, it is possible to represent it uniquely as an element of the Special Linear group SL(3) by choosing ζ such that det(H) = 1, namely

H = H det(H) -1 3 ∈ SL(3).
For the remainder of this paper, we assume that the calibrated homograpy is appropriately scaled in such a way H ∈ SL(3).

D. Homography Kinematics and Measurements

Denote the angular velocity and linear velocity of the target frame {T } with respect to the inertial frame {I} expressed in the target frame {T } by Ω T and V T . The target kinematics of (R T , ξ T ) are given by:

ṘT = R T Ω T × , ξT = R T V T . ( 10
)
Let Ω denote the angular velocity of the body-fixed frame expressed in {B} and let V denote the linear velocity of the body-fixed frame {B} with respect to the inertial frame {I} expressed in the frame {B}. The kinematics of (R, ξ) are given by:

Ṙ = RΩ × , ξ = RV. ( 11
)
It follows that the kinematics of the relative attitude and position of the body-fixed frame {B} with respect to the target frame {T } expressed in {T } can be written as

Ṙ = R Ω× , ξ = R V (12)
where

Ω× = Ω × -( R Ω T ) × , (13a) V = R Ω T × ξ + V -R V T . ( 13b
)
Lemma 1. Consider a camera moving with kinematics (11) viewing a moving planar scene whose kinematics are defined by [START_REF] Chen | Homographybased visual servo tracking control of a wheeled mobile robot[END_REF]. Let H : {B} → {T } denote the Euclidean homography matrix. The group velocity Ū ∈ sl(3) induced by the relative motion between the camera and the target plane is such that

Ḣ = H Ū , with Ū = Ω× + V η d - η V 3d I. (14)
Proof: See Appendix A. Note that the group velocity Ū induced by the relative motion between the camera and the target plane depends on the normal to the surface η and the orthogonal distance d, that define the scene geometry at time t, as well as the angular velocity of the target plane Ω T . The scene parameters and the velocity of the target are usually unmeasurable and they have to be taken into account during the observer design. To this purpose we rewrite the group velocity as

Ū = Ω × + Γ ∧ , ( 15
)
where Ω is typically the measurable part, which can be obtained from the set of embedded gyros and can be directly used in the observer as feed-forward term. While the remainder of the group velocity Γ ∧ represents the unmeasurable part that has to be estimated.

In order to cope with the unmeasurable group velocity, we assume that Γ ∧ can be immersed into a finite-dimensional and observable linear system ( [START_REF] Byrnes | Structurally stable output regulation of nonlinear systems[END_REF], [START_REF] Khalil | Robust servomechanism output feedback controllers for feedback linearizable systems[END_REF]) in such a way it can be written in the form

Γ ∧ = (Cγ) ∧ , γ = Sγ ( 16
)
where the known matrix C ∈ R 8×m is a full row rank matrix with m ≥ 8, γ ∈ R m and the known matrix S ∈ R m×m is skew symmetric, namely S = -S .

Remark. The eigenvalues of a real skew-symmetric matrix are all purely imaginary and they come in conjugate pairs. In case of odd dimension there is an additional unpaired null eigenvalue. The system in ( 16) is stable in the Lyapunov sense. Indeed considering the Lyapunov candidate L = γ 2 , and differentiating it along the solutions of (16) one obtains L = 2γ Sγ = 0. It follows that any autonomous observable system of the form

ẋ = Ax, A ∈ R m×m , x ∈ R m y = C y x, C y ∈ R 8×m , y ∈ R 8
is equivalent to (16) iff all eigenvalues of A have zero real part and multiplicity one in the minimal polynomial. The considered assumption on Γ ∧ is of course not valid for all types of relative motions between the camera and the target plane. It, however, allows us to cover many possible camera/target trajectories encountered in practice such as (near) circular or sinusoidal motions.

We consider a set of n point measurements p i ∈ P 2 associated to a group action of SL(3) onto the projective space

P 2 p i = H -1 pi |H -1 pi | , i = {1, . . . , n}, (17)
representing calibrated image points projected onto the unit sphere.

III. NON-LINEAR OBSERVER DESIGN ON SL(3)

The goal of the non-linear filter is to provide an estimate Ĥ ∈ SL(3) given a collection of n measurements p i to drive the error H := ĤH -1 to the identity element of the group. To this purpose, we define the estimates pi ∈ P 2 of p i as

pi := Ĥ-1 pi Ĥ-1 pi , i = {1, . . . , n}. (18)
The estimates e i of pi then is defined as

e i := Ĥp i Ĥp i = Hp i Hp i , i = {1, . . . , n}. (19)
Definition 1. A set M n of n(≥ 4) vector directions pi ∈ P 2 is called consistent, if it contains a subset M 4 ⊂ M n of 4 constant vector directions such that all its vector triplets are linearly independent.

Theorem 1. Let H denote the Euclidean homography (9) and consider the kinematic system in (14) along with (15) and ( 16).

Consider the following nonlinear filter

Ḣ = Ĥ(Ω × + Γ∧ ) + k p ∆ Ĥ (20a) Γ∧ = (C γ) ∧ (20b) γ = Sγ + k I C Q sl3 n i=1 ϕ i (20c)
with k p and k I some positive gains and

∆ := n i=1
π ei pi e i ∈ sl(3)

ϕ i := vec ∨ P sl3 ( Ĥ pi )( Ĥ-1 e i ) +
-e i pi vec ∨ P sl3 ( Ĥ e i )( Ĥ-1 e i ) .

(21) Assume that the measured angular velocity Ω of the planar target, the group velocity Γ (equivalently γ) in ( 16) and homography matrix H are bounded.

Then, if the set M n of measured directions pi is consistent, the equilibrium ( H, γ) = (I, 0), with H := ĤH -1 , γ := γ -γ, of the error system is locally exponentially stable.

The proof of Theorem 1 is given in Appendix B. The positive gains k p and k I in (20) share the same effects of a proportional gain and an integral gain for a linear system, respectively. In particular k p directly relates the dynamics of Ĥ to the innovation, while the gain k I represents the integral term of the innovation projected into R 8 in the inertial frame.

Remark. The constant image points pi represent known reference points which can be extracted from a stored reference image of the planar scene or extracted directly during the initialization process of the observer, namely the first image taken by the camera. In the latter case one has pi = p i (0), ∀i = {1, . . . , n}, whereas e i represent the estimates of the reference points pi based on the current measurements p i . Roughly speaking, in a classical linear Luenberger observer the error i (p i -e i ) would represent the output prediction error (innovation).

IV. SIMULATION RESULT

In this section we illustrate the performance of the proposed observer through simulations. The camera is assumed to be a stationary camera. The target plane trajectory, depicted in Figure 3, is chosen as a Lemniscate of Gerono with zero angular velocity and linear velocity given by:

V T =   -12π sin(6πt) 12π cos(12πt) 24π cos(12πt)   .
Note that the group velocity Γ ∧ in [START_REF] Francis | The internal model principle of control theory[END_REF] induced by the Gerono trajectory cannot be immersed into a finite linear system. However, due to the fact that Γ ∧ is a periodic function it is possible to decompose it into Fourier series (see Figure 4).

The initial homography H is initialized as while the initial estimated homography Ĥ(0) is set to the identity element of the group. The filter parameters are k p = 50 and k I = 10.

H(0) =   2
Considering the spectrum of the group velocity induced by the motion of the planar target, depicted in Figure 4, we implemented in the internal model unit the first four harmonics for the basis B 1 , B 5 , B 7 , B 8 defined in [START_REF] Astolfi | A high-gain nonlinear observer with limited gain power[END_REF].

For comparison purpose a second simulation was run considering the observer proposed in [START_REF] Hamel | Homography estimation on the special linear group based on direct point correspondence[END_REF], using the same value k p = 50 for the proportional term. The results reported in Figures 56show respectively the time behavior of the Frobenious norm I -H F for the observer presented in [START_REF] Hamel | Homography estimation on the special linear group based on direct point correspondence[END_REF] and the proposed observer (termed OSC in order to highlight the oscillatory nature of the terms in the Lie algebra), the Frobenious norm Γ ∧ -Γ∧ F , and the group velocity vs. the estimated group velocity. Figure 5 clearly shows that the proposed observer outperform the observer in [START_REF] Hamel | Homography estimation on the special linear group based on direct point correspondence[END_REF]. Plots show that, in steady state, the velocity Γ∧ of the observer practically converges to Γ ∧ and the estimated homography Ĥ practically converges to the real homography H. In red the proportional observer (P) proposed in [START_REF] Hamel | Homography estimation on the special linear group based on direct point correspondence[END_REF] and in blue the observer (OSC) proposed in this paper.

In Figure 7 it is shown the state of the internal model unit relative to the basis B 8 in order to track the desired group velocity Γ ∧ .

V. EXPERIMENTAL RESULT A. Experimental implementation aspects

The proposed non-linear observer has been implemented in C++ with OpenCV library. Due to real-time constraint, the time consuming image processing part of the observer has been implemented taking advantages of modern GPU parallel computing. In particular, in view of embedding the filter on small-scale unmanned vehicles on the Jetson NVIDIA family, CUDA API have been used. The algorithm, whose flowchart is depicted in Figure 8, can be summarized as follows. The process is initialized detecting feature and extracting the descriptor of a reference image (red boxes Figure 8). As soon as a new image arrives, it is transformed with a perspective transformation (OpenCV's cuda::warpPerspective) based on the homography estimate. The algorithm detects key-points on the current warped image and extracts the features' descriptor. The descriptor of the current image is matched with the descriptor of the reference image (Figure 9). Those matches, however, may contain many outliers which are removed with the outliers algorithm presented in [START_REF] Hua | Featurebased Recursive Observer Design for Homography Estimation and its Application to Image Stabilization[END_REF]. The remaining inliers, that lie on the image plane, are projected onto the unit sphere [START_REF] Gauthier | Deterministic Observation Theory and Applications[END_REF]. From these measurements the homography estimate is updated by iterating the observer equations 250 times per video frame. In this correction step, special attention has been paid on the integration of the observer equations in [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. In particular a leapfrog integrator has been implemented for the numerical integration of the oscillators bank in the lie algebra sl(3) of the observer. Note that feature detection and description are essential components for the algorithm. Several feature detectors, such as scale invariant feature transform (SIFT) [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], features from accelerated segment test (FAST) [START_REF] Rosten | Fusing points and lines for high performance tracking[END_REF], oriented FAST and rotated BRIEF (ORB) [START_REF] Rublee | Orb: An efficient alternative to SIFT or SURF[END_REF], speeded-up robust features (SURF) [START_REF] Bay | Speeded-up robust features (SURF)[END_REF] and descriptors such as ORB, SIFT and SURF have been proposed in literature. Choosing the detector-descriptor combination depends on the particular application and is usually a trade-off between accuracy, repeatability, robustness and speed. Since the major constraint for this application is the time it has been decided to use the combination FAST-ORB.

B. Experimental setup

In order to experimentally validate the proposed observer we present an application for the homography estimate of an oscillating cardboard target in challenging conditions such as partial target occlusions, changing in light condition and specular reflections. The platform used for the experiment, depicted in Figure 10, is made of a cardboard box whose upper flap is actuated by a small servo motor via an Arduino Mega 2560. The camera used is a Basler acA 1200-200uc featuring a 6.1 mm × 4.9 mm sensor, maximum resolution 1280 px × 1024 px and maximum frame rate of 203 fps. An Optitrack motion capture system is used together with markers mounted to the camera and to the cardboard in order to provide homography ground truth measurements starting from the knowledge of the full pose of the camera and the planar target.

C. Experimental results

The experiment has been performed using the Basler camera which is held stationary recording at 50 fps while looking at the oscillating target. In an ideal scenario, the target platform should oscillate with a frequency of 0.83 Hz only on the pitch angle. However, due to the simple design of the handmade platform, the target oscillates along all the three axes, as depicted in Figure 11. Since the platform is tracking the desired sinusoidal trajectory by means of position feedback only, and due to the fact that the center of mass of the target is located slightly below the rotation axis, the rising and the falling edge of the sinusoidal trajectory are not symmetrical. As a consequence, the motion of the target is still periodical but with a much richer frequency spectrum (see Figure 11, cyan box).

We proceed by presenting a comparison of the homography observer proposed in [START_REF] Hamel | Homography estimation on the special linear group based on direct point correspondence[END_REF] with respect to the filter formulated in this paper. In order to have a fair comparison between the algorithms, the proportional term of the three observers is chosen as k p = 80, the integral term of the proportionalintegral filter is chosen equal to the oscillator gain of the internal model-based observer and set to k I = 80. In the internal model-based observer only the nominal oscillation frequency f n = 0.83 Hz of the target has been implemented in the bank of oscillators defined on the Lie algebra sl(3). In particular denoting by

S 1 := 0 -2πf n 2πf n 0 , C 1 := 0 1 ,
the matrices of the implemented observer have the following form

S = diag(S 1 , S 1 , S 1 , S 1 , S 1 , S 1 , S 1 , S 1 ) ∈ R 16×16 , C = diag(C 1 , C 1 , C 1 , C 1 , C 1 , C 1 , C 1 , C 1 ) ∈ R 8×16 .
Figures [START_REF] Bastin | Stable adaptive observers for nonlinear time-varying systems[END_REF][START_REF] Bay | Speeded-up robust features (SURF)[END_REF][START_REF] Benhimane | Homography-based 2D visual tracking and servoing[END_REF][START_REF] Bin | Robust internal model design by nonlinear regression via low-power high-gain observers[END_REF][START_REF] Byrnes | Output regulation of Uncertain Nonlinear Systems[END_REF][START_REF] Byrnes | Structurally stable output regulation of nonlinear systems[END_REF][START_REF] Chaumette | Visual servo control part I: Basic approaches[END_REF][START_REF] Chaumette | Visual servo control, part II: Advanced approaches[END_REF][START_REF] Chen | Homographybased visual servo tracking control of a wheeled mobile robot[END_REF][START_REF] Davison | The robust control of a servomechanism problem for linear time-invariant multivariable systems[END_REF][START_REF] Marco | Homography estimation of a moving planar scene from direct point correspondence[END_REF] show the comparison between the three algorithms: namely proportional (P), proportional-integral (PI), and internal model-based (OSC). To this regard, we performed an experiment where:

• From 0-25 seconds (Figure 2 frames 38-293) the conditions are ideal, except for the presence of specular reflections. Figure 12 shows the time behavior of the Frobenious norm I -H F for the three implemented algorithms. The figure clearly shows better performance of the internal model based observer with respect to [START_REF] Hamel | Homography estimation on the special linear group based on direct point correspondence[END_REF] in all the conditions except when the assumption on the periodicity of the target trajectory is strongly violated. APPENDIX A PROOF OF LEMMA 1 Considering (9), the time derivative of the homography is given by

Ḣ = ζ Ṙ + ξη + ξ η d - ḋ ξη d 2 + ζ ζ H.
Recalling [START_REF] Byrnes | Structurally stable output regulation of nonlinear systems[END_REF], [START_REF] Chaumette | Visual servo control part I: Basic approaches[END_REF] it is straightforward to verify that ḋ = -η V , η = η Ω× and bearing in mind (13a), (13b) one has

Ḣ= ζ R Ω× + R V η + ξη Ω× d + η V ξη d 2 + ζ ζ H = ζ R + ξη d Ω× + ζ R + ξη d V η d + ζ ζ H = H Ω× + V η d + ζ ζ I . (22

APPENDIX B PROOF OF THEOREM 1

We recall that the matrix representation of the composition of the linear maps (vect •∧) is denoted by the full columns rank matrix J∧K ∈ R 9×8 . Let Jvec ∨ K ∈ R 8×9 denotes the matrix representation of the inverse of the map (vect •∧), namely

vect(v ∧ ) = (vect •∧)(v) = J∧Kv, vec ∨ (v ∧ ) = (vect •∧) -1 (vect(v ∧ )) = Jvec ∨ K vect(v ∧ ), (23)
for any v ∈ R 8 . The operator ⊗ denotes the usual Kroneker product. Let us introduce the following properties which are instrumental for the proof of the main theorem of this paper.

A. Property 1

For any H ∈ SL(3) and v ∈ R 8 , the matrix representation of the composition of the linear maps (vec

∨ • Ad H ) vec ∨ (Ad H (v ∧ )) := JAd H Kv is given by JAd H K = Jvec ∨ K(H -⊗ H)J∧K. (24)
Proof. From the definition of J∧K and Jvec ∨ K in (23) one has

vec ∨ (Hv ∧ H -1 ) = Jvec ∨ K vect(Hv ∧ H -1 ) = Jvec ∨ K(H -⊗ H) vect(v ∧ ) = Jvec ∨ K(H -⊗ H)J∧Kv = JAd H Kv B. Property 2 For any A ∈ R 3×3 , B ∈ R 3×3 , one verifies J∧KJvec ∨ K(A ⊗ B)J∧K vec ∨ (U ) = = J∧KJvec ∨ K(A ⊗ B) vect(U ) = J∧KJvec ∨ K vect(BU A ) = vect(BU A ) = (A ⊗ B) vect(U ) = (A ⊗ B)J∧K vec ∨ (U ) which implies J∧KJvec ∨ K(A ⊗ B)J∧K = (A ⊗ B)J∧K. (25)

C. Theorem 1 Proof

By considering the Euclidean homography ( 14), with kinematics ( 15), ( 16) along with the observer [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF], the dynamics of the error system are given by

Ḣ = -[Ad Ĥ (C γ) ∧ ] H + k p ∆ H γ = Sγ -k I C Q sl3 n i=1 vec ∨ (Ad Ĥ ∆). (26)
To prove that the origin of the error system ( H, γ) = (I, 0) is locally exponentially stable, it suffices to show that the origin of the linearized error system is exponentially stable. The proof is based on Theorem 1 in [START_REF] Loria | Uniform exponential stability of linear timevarying system: revisited[END_REF] which establishes sufficient conditions for the uniform exponential stability of the origin of a linear time-varying system having the following standard form ẋ θ = A(t) B(t) -C(t) 0

x θ . (27)
In order to write the linearized system in the standard form [START_REF] Loria | Uniform exponential stability of linear timevarying system: revisited[END_REF] , which implies x Px ≥ 0. We show that P is positive definite by contradiction. Assume that there exists a non null vector x such that x = 0 ∈ ker P,

which in turn implies and recalling that π pi is a projector it yields x∧ pi = λ i pi .

Since the measurement set is consistent, one can consider (without loss of generality) that (p 1 , p2 , p3 ) are three non collinear eigenvectors of x∧ associated with the eigenvalues λ i for i = {1, 2, 3}. Moreover due to the consistency of the set it can be shown that there exists a constant direction pk from the set {p 4 , . . . , pn } such that: which implies x Qx ≥ 0, and using the fact that the set M n is consistent one can ensure that Q is positive definite. Finally, it is straightforward to verify that B is a full rank matrix, which implies that B B is positive definite which in turn implies that the term B is persistent exciting. Therefore, we conclude that all the conditions of Theorem 1 in [START_REF] Loria | Uniform exponential stability of linear timevarying system: revisited[END_REF] are satisfied and thus the set S is locally exponentially stable. This in turn completes the proof.

Fig. 1 .

 1 Fig. 1. Representation of a point of the planar surface in the inertial frame {I}, body-fixed frame {B} and target frame {T }.

Fig. 2 .

 2 Fig. 2. The red, blue and green polygon represents the region of interest tracked by the proportional, proportional-integral and internal model-based filter respectively.

Fig. 3 .

 3 Fig. 3. Target Plane Trajectory.

Fig. 4 .

 4 Fig. 4. Frequency spectrum of the group velocity Γ∧.

Fig. 5 .

 5 Fig. 5. Time behaviour of the Frobenious norm I -H F and Γ∧ -Γ∧ F .In red the proportional observer (P) proposed in[START_REF] Hamel | Homography estimation on the special linear group based on direct point correspondence[END_REF] and in blue the observer (OSC) proposed in this paper.

Fig. 6 .

 6 Fig. 6. Estimated group velocity Γ∧ (blue lines) and true group velocity Γ∧ (red lines). Note that each plot represent an element of the group velocity matrix.

Fig. 7 .

 7 Fig. 7. Bank of oscillators relative to the eighth generator of the Lie algebra sl(3). In blue the first four harmonics (6, 12, 18, 24 Hz) implemented in the observer and in cyan the resultant.

Fig. 8 .

 8 Fig. 8. Flowchart of the implemented algorithm.

Fig. 9 .

 9 Fig. 9. Key-points detection, description and matching.

Fig. 10 .

 10 Fig. 10. Oscillating Target Platform.

Fig. 11 . 4 PIFig. 12 .

 11412 Fig. 11. Roll, pitch and yaw Euler angles of the oscillating platform.

  we consider the following time-varying change of coor-dinates γ = Qγ, Q = -QS, with Q ∈ SO(m), yielding Ḣ = -[Ad Ĥ (C Q γ ) ∧ ] H + k p ∆ H, γ = -k I QC Q sl3 n i=1 vec ∨ (Ad Ĥ ∆).(28)Let us define x ∧ , with x ∈ R 8 , and θ ∈ R m the first order approximation of H, γ around the equilibrium point (I, 0)H (I + x ∧ ) = (I + 8 k=1 x k B k ), γ = θ. ⊗ π pi ) (p i ⊗ π pi )J∧Kx = vect(x ∧ ) n i=1 (p i ⊗ π pi ) (p i ⊗ π pi ) vect(x ∧ ) = n i=1 p i x ∧ π pi x ∧ pi = n i=1 e i -pi2 H=I

π pi x∧ pi = 0 .

 0 It follows that for each i = {1, . . . , n} one has x∧ pi ∈ ker(π pi )

pk = ẙk |ẙ k | , where ẙk = 3 i=1α 3 i=1 α i pi = 3 i=1 λ i α i pi = 3 i=1λ

 3333 i pi , α i ∈ R/{0}, i = {1, 2, 3} which implies that x∧ k α i pi .From the equation above it is straightforward to verify thatλ 1 = λ 2 = λ 3 = λ k . Since x∧ ∈ sl(3) one has tr(x ∧ ) = 3 i=1 λ i = 0,which along with the equation above implies λ 1 = λ 2 = λ 3 = λ k = 0, which in turn implies x = 0 and this contradicts x = 0 in[START_REF] Nguyen | Inertial-aided homography-based visual servo control of autonomous underwater vehicles without linear velocity measurements[END_REF]. It follows that the kernel of the matrix P is trivial and P is positive definite.For the matrix Q, one verifiesx Qx 2k I k p = x J∧K n i=1 (p i p i ⊗ π pi ) n i=1 (p i p i ⊗ π pi )J∧Kx = vect n i=1 (π pi x ∧ pi p i ) vect n i=1(π pi x ∧ pi p i ) = tr(∆ ∆) H=I ,

ACKNOWLEDGMENT

This research was supported by the French DGA-Rapid project "Alcyon", the French EQUIPEX Robotex project, the French Astrid CONGRE project ANR-18-ASTR-0006 and the Australian Research Council through the "Australian Centre of Excellence for Robotic Vision" CE140100016.

The authors also thank C. Samson (INRIA and I3S) for his constructive feedback and advice on this work.

A first-order approximation of the output errors e i given by [START_REF] Hamel | Homography estimation on the special linear group based on direct point correspondence[END_REF], considering the equation above, can be written as

and by neglecting high order terms one has e i pi + x ∧ pi -pi p i x ∧ pi .

Substituting the expression of the output errors e i [START_REF] Malis | 2-1/2D visual servoing[END_REF] in the innovation term ∆ one obtains

(31) Thus, the first-order error dynamics for the error system are given by

Bearing in mind Property 1 in [START_REF] Isidori | Output regulation for nonlinear systems[END_REF], it yields

Proceeding in a similar way, the first order approximation of γ is given by

and recalling the definition of Q sl3 and bearing in mind Property 2 in (25) one obtains

It follows that the linearized system is in standard form [START_REF] Loria | Uniform exponential stability of linear timevarying system: revisited[END_REF] with

It is straightforward to verify that the matrices

are symmetric and they satisfy the required relations PB = C and Q = -A P -PA of Theorem 1 in [START_REF] Loria | Uniform exponential stability of linear timevarying system: revisited[END_REF]. It remains to prove that the matrices P and Q are positive definite, to this purpose note that the matrix P is associated to the quadratic form e i -pi 2 at the origin, indeed one has