Bassem Dahroug

Jean-Antoine Séon

Brahim Tamadazte

Nicolas Andreff

Ali Oulmas

Tiantian Xu

Stéphane Régnier

Some examples of path following in microrobotics

Keywords: microrobot control, biomedical robotics, path following, visual servoing, laser steering, magnetic manipulation

After reviewing the various motion generation strategies (waypoint sequence, trajectory tracking and path following), this paper suggests that path following is probably the best strategy at small scales. This is illustrated by three examples coming from biomedical applications: endoscopic laser steering, magnetic manipulation of microswimmers and navigation of a surgical tool under anatomical constraints.

I. INTRODUCTION

Small scale creates specific challenges to robotics related to a different ratio between volumic forces (namely, gravity and dynamics) and surface forces (namely, friction) than at conventional scale [START_REF] Chaillet | Microrobotics for Micromanipulation[END_REF]. One of them is related to motion control. Indeed, at small scale, inertia is often neglectable, which has two effects. The first one is that small objects can reach a very high velocity almost immediately, the second one is that friction (stiction at the solid/air interface or fluidic drag) becomes a predominant disturbance. Therefore, the conventional motion control strategies might have to be reconsidered.

This issue was encountered when tuning the controller for the MagPier, on the occasion of the NIST mobility challenge at ICRA 2011 [START_REF] Ivan | Nist and ieee challenge for magpier: The fastest mobile microrobots in the world[END_REF]. Indeed, the mobility challenge consisted in a slalom through 4 gates, and the MagPier was a magnetic particle driven by a tunable magnetic field. Therefore, a "trajectory" (see below for a proper definition) passing through the gates was defined and conventional robot control was applied, but it would not work.

What happened ? In fact, the trajectory was defined as a sequence of waypoints and designed the controller as a point to point step response, with a simple PID control providing the force towards the next waypoint, this force being then converted into currents into the magnetic coils by appropriate transduction. Thereby, the MagPier would stay stuck on the substrate until the integral term would break the stiction, then flew until the distance to the closest waypoint shrunk so small that the next waypoint would be set as a target. At least, this was the expected behaviour. However, when breaking stiction, the magnetic force eventually was so large that the MagPier flew over the waypoint over to a further one. The controller being simple, not to say stupid, was therefore unable to stabilize the MagPier on the desired waypoint before switching to the next one. So, an ever more stupid motion controller was designed: "fly towards the next waypoint, update your position, determine which is the closest waypoint ahead on the sequence, and start again". And this stupid algorithm worked so well that the MagPier won the challenge twice in a row.

What did we do ? Actually, we rediscovered path following, a motion control strategy that is time independent and is essentially used for automated guided vehicles. It somehow differs with the motion strategies used in industrial robotics, where controller are defined with a dependence on time, as it will be explained in Section II. Turning later to biomedical applications, we found that path following could be profitably be used. Indeed, it allows to decouple laser tissue/interaction from anatomical constraints in laser surgery (Section III), it allows to derive a proper controller for microswimmers subject to non-holonomic constraints (Section IV) and it can even be coupled with remote center of motion (RCM) constraints in keyhole surgery (Section V).

The contribution of this paper is to trace our experience in path following at small scale through these examples, which explains the journalistic style used throughout the paper and the egocentrism of the references.

II. MOTION CONTROL STRATEGIES

A. Waypoint sequence

The first motion control strategy consists in defining a sequence of intermediate positions (waypoints) to be reached by the robot as successive step responses:

S = {P * k ∈ ℜ 3 |∀k ∈ {1, ..N }}.
More complex motion is achieved with a sequence over SE(3) instead of ℜ 3 to impose constraints on the orientation. It is thus extremely simple for the user. It is also very easily programmable since proportional control u(t) = K p * (P(t) -

P * k ) (1)
is often enough in the absence of strong disturbances if the system is linear or close to be linear. When disturbances occur, this controller is most often replaced by a PID controller, another linear controller. However, since the microrobot dynamics are often non-linear, tuning of the linear gains for any step might push those gains downwards. More important, this strategy provides a velocity profile, with exponentially decaying velocities from one waypoint to the other one. This profile is not really smooth, since the robots would alternately start and stop. This might even be troublesome for mobile microrobots in presence of friction.

To overcome these drawbacks, two classes of techniques can be found: interpolation between waypoints under velocity constraints (that gives rise to trajectory tracking) and appropriate accountance of differential geometry (that gives rise to path following).

B. Trajectory tracking

Trajectory tracking is the major motion control strategy in industrial robotics. A trajectory is defined by a set of N (time, position) couples:

{(t k , P * (t k )) ∈ ℜ × ℜ 3 |∀k ∈ {1, ..N }}.
Again, more complex trajectories can be defined over SE(3) instead of ℜ 3 to include orientation constraints. In general, a few number of such couples are defined, and then interpolation of various orders fills in the gap to densify the set: T = {P * (t) ∈ ℜ 3 |∀t ∈ ℜ} while preserving some smoothness property (usually C 1 or C 2 ) of the function P * (t) with respect to time.

The advantages of this motion control strategy are many. Among others, it may provide optimality in the overall traversal time, it may take into account the robot dynamics and it is rather simple to implement with a proportional+feedforward controller:

u(t) = K p * (P(t) -P * (t)) + dP * (t) dt (2)
Again, this controller can be extended to a PID+feedforward to account for small non-linearities and disturbances. The drawbacks of this strategy are that it requires an expert roboticist for choosing the time constraints, whereas the nonroboticist user generally does not have the slightest hint on the latter. Moreover, specifically to low Reynolds swimming, where the swimming velocity varies with the fluid viscosity, such a control strategy is not robust to time delays and may thus create geometric deviations.

C. Path following

A path is defined as a smooth curve, at least C 1 , with respect to the distance s ∈ ℜ traveled along the curve:

Γ = {P(s) ∈ ℜ 3 |∀s ∈ ℜ}. It is thus independent from time.
And so should the controller for its following be. Therefore, this motion control strategy decouples geometry from velocity, as a car can drive the same path on a motorway at various velocities.

The main drawback of path following is that it seems rather complex to implement. However, the advantages of it are many. From a technical viewpoint and by definition, path following is robust to time delays along the path and fully respects the desired geometric path. More important, it remains intuitive for the user to define the path (provided that the human-robot interface is adequately designed). From an application viewpoint, decoupling the geometric path from the velocity is particularly useful to biomedical applications of microrobotics. Indeed, for biological applications, respecting the geometric path is crucial to prevent contamination. For clinical applications, path-velocity decoupling garantees the respect of the surgical planning and allows the surgeon to share the control with the robot, the robot bringing accuracy and the surgeon keeping his hands on the tool progression.

III. ENDOSCOPIC LASER STEERING

The first example in this paper arised in the FP7 µRALP project, dedicated to microrobot-assisted laser phonosurgery. In this context, the clinical specifications are that the surgeon would use the endoscopic image to delineate the tumor to be removed from the vocal folds with a laser and that the laser would move at a constant speed. Indeed, moving too slow, the laser would burn the tissues instead of incising them; moving too fast, it would simply have no effect.

Therefore, we decided to take advantage of the path-velocity decoupling within path following. Since the laser was observable in the endoscopic image, in which the desired path is defined by the surgeon, the problem was to drive a 2D point along a 2D path. A similar problem was already solved in the literature for automated guided vehicles, where an additional non-holonomic constraint is at stake: the vehicle can not move sideways. In order to have smooth motion of the laser (and to avoid zig-zaging scars), we deliberately applied the same nonholonomic constraint to the laser, although it was steered by a mirror and was physically not subject to such a constraint [START_REF] Séon | Decoupling path following and velocity profile in vision-guided laser steering[END_REF].

Thereby, the state of the laser spot is defined by (s, d, θ e ), where s is the curvilinear abscissa of the closest point h p to the laser spot on the path, d is the shortest distance from the laser spot to the path and θ e is the angle between the laser spot velocity u 1 and the vector x s , tangent to the path in h p (Figure 1). The virtual kinematics of the laser spot were then given by the unicycle kinematics: where u 1 = u 1 is the norm of the laser spot velocity, u 2 = θe is the norm of the angular velocity and C(s) is the curvature of the path in h p . Note that u 1 and u 2 are the two control inputs, u 1 being dedicated to the laser spot velocity and u 2 being in charge to servo the laser spot on the path at the velocity given by u 1 .

ṡ = u 1 1 -dC(s) cos θ e (3)
ḋ = u 1 sin θ e (4) θe = u 2 -ṡC(s) (5)
In that scheme, the non-roboticist user only needs to define the path Γ and the laser velocity u 1 . For instance, the former can be drawn with a stylet on a touchpad (Figure 2), while the latter can be intuitively controlled by foot pedals as in a car. It is the role of the developer to compute u 2 in function of the detected position of the laser spot in the image and those user inputs. The solution to this problem is elegantly given in [START_REF] Morin | Motion control of wheeled mobile robots[END_REF] by an appropriate change of coordinates that casts the kinematics in a 3-state/2-inputs cascade system:

ż1 = ṡ (6) ż2 = ṡz 3 (7) ż3 = v 2 (8)
where ṡ depends on the state and u 1 while u 2 can be deduced from the definition of v 2 in the change of coordinate:

v 2 = (-ḋC(s) -d ∂C(s) ∂s ṡ) tan θ e +(1-dC(s))(1 + tan 2 θ e )u 2 (
9) and from the stable control established for such a cascade system:

v 2 = -v 1 γ 2 z 2 -|v 1 | γ 3 z 3 ( 10
)
where γ 2 and γ 3 are positive gains. This control was associated to a parallel pan-tilt microrobot to move the laser and tested both on a testbed and in a human cadaver trial [START_REF] Renevier | Endoscopic laser surgery: Design, modeling and control[END_REF] with an accuracy of 30µm on the testbed and of 84µm in the preclinical test (Figure 3).

IV. MAGNETIC MANIPULATION OF A MICROSWIMMER

Microswimmers, be they helical [START_REF] Tottori | Magnetic helical micromachines : Fabrication, controlled swimming, and cargo transport[END_REF] or flagellar [START_REF] Khalil | Magnetic-based motion control of sperm-shaped microrobots using weak oscillating magnetic fields[END_REF], [START_REF] Ye | Rotating magnetic miniature swimming robots with multiple flexible flagella[END_REF], convert a time-varying uniform field field into motion thanks to their tail. This creates a non-holonomic constraint since the Fig. 3. Pre-clinical trial of endoscopic laser steering [START_REF] Renevier | Endoscopic laser surgery: Design, modeling and control[END_REF] microswimmer can not move sideways. Moreover, the velocity amplitude of a helical swimmer can not be reliably achieved in many cases. The first reason for that is the well-known cutoff frequency of the velocity amplitude [START_REF] Honda | Micro swimming mechanisms propelled by external magnetic fields[END_REF], which imposes a strong saturation of the available ampli-tude. As we have shown in [START_REF] Xu | Modeling and swimming property characterizations of scaled-up helical microswimmers[END_REF], this cut-off frequency depends on many uncontrolled factors (microfabrication quality and environmental viscosity, for the main ones). The second reason is related to the potential stiction of the helix onto the substrate, which creates a distortion in velocity amplitude. The third reason is related to the viscosity at low scales (through the Reynolds number). Indeed, at low scales, we are not sure that the fluid is really homogeneous (inherently or due to the presence of micro or nano-dirt in a non totally clean room). This is especially true when dealing with biological fluids (blood, cell culture) where the medium is definitely heterogeneous or, at least, where the cells represent mobile obstacles that need to be pushed aside to allow for the helical swimmer to follow the expected trajectory.

For the non-roboticist user, the task for a microswimmer is usually to cargo a cell from one point to another one. However, with appropriate human-machine interface and with advanced control, the user will become able to define more complex tasks, such as navigating in a Petri dish through a complex map defined by several cell clusters separated by empty spaces. The non-holonomic constraint makes it natural to apply the controller defined in the previous section. This allowed us to accurately guide a helical swimmer in closed loop in a Fig. 4. Open-loop vs. closed loop control in planar path following [START_REF]Planar path following of 3d steering scaled-up helical microswimmers[END_REF] Fig. 5. Helical swimmer modelling [START_REF] Oulmas | Chained formulation of 3D path following for nonholonomic autonomous robots in a Serret-Frenet frame[END_REF] horizontal plane (Figure 4), while the height of the swimmer was precompensated in open loop [START_REF]Planar path following of 3d steering scaled-up helical microswimmers[END_REF].

The planar controller was then extended to the 3D case, by extending the unicycle kinematics over a 5 degrees-offreedom state defined as (s, d y , d z , θ de , θ ie ) where s is the curvilinear abscissa of h p (the orthogonal projection of the swimmer position onto the path), d y and d z are the projection of the distance from the swimmer to the path on two axes perpendicular to the path tangent vector, and finally θ de and θ ie are the two angles defining the angular error between the swimmer velocity and the path tangent vector (Figure 5). Then, an appropriate change of coordinates allowed us to transform the swimmer model into a 5-state/3-inputs cascade system, from which we were able to derive a controller and prove the stability of the latter [START_REF] Oulmas | Chained formulation of 3D path following for nonholonomic autonomous robots in a Serret-Frenet frame[END_REF]. With that, we were able to follow a helix (Figure 6) with an accuracy better than 2% of the swimmer length [START_REF]Closed-loop 3D path following of scaled-up helical microswimmers[END_REF]. Changing only the magnetic actuation mode but not the 3D path following controller [START_REF]3D closed-loop motion control of swimmer with flexible flagella at low Reynolds numbers[END_REF], we were also able to navigate a swimmer with a flexible flagella along an arbitrary path (Figure 7).

V. NAVIGATION OF A SURGICAL TOOL UNDER ANATOMICAL CONSTRAINTS

The last example comes from another surgical context, this time in the case of middle-ear surgery. Again, the clinical specifications request accurate 3D path following for ablating Fig. 6. Following a helix with a helical swimmer [START_REF] Oulmas | Chained formulation of 3D path following for nonholonomic autonomous robots in a Serret-Frenet frame[END_REF] Fig. 7. Following an arbitrary path with a flexible flagella [START_REF]3D closed-loop motion control of swimmer with flexible flagella at low Reynolds numbers[END_REF] pathologic tissues around the ossicles. In addition, as in many other surgical case, the access to the middle ear is restricted to a single point while realizing the path following (Figure 8).

The definition of the 3D path can be done from any 3D surgical planner. However, the specifications do not mention any constraint on the velocity. In fact, the exact time to finish the path is not a factor of significant nature for the surgeon. What really matters is that the instrument follows exactly the intended path. The surgeon may also need to change the tool velocity independently from the path shape, size or curvature. The instrument velocity is indeed dependent on the tissue-tool interaction which could be variable due to the potential inhomogeneity of the tissue. In addition, the tool Fig. 8. Path following under single point access constraint [START_REF] Dahroug | Visual servoing controller for time-invariant 3D path following with remote centre of motion constraint[END_REF] velocity profile may change when the surgeon scans a ROI, or moves from one ROI to another for inspection. For the above reasons, the path following controller is a good option for our application; since it separates the geometric curve from the velocity profile.

Note also that this is coherent with surgical planners. Indeed, as far as we know, surgical planners restrict their work to path definition and do not time stamp the path into a trajectory.

On the contrary to the first two examples, non-holonomic constraints were neither physically imposed nor virtually chosen. Only the position of the tool tip was considered in the control. Thereby, the state of the system is given by (s, d) where d = ph p is the vector, not scalar, distance from the path to the tool tip. After geometric analysis in 3D, the time derivative of the lateral error becomes;

ḋ = I 3×3 - x s x T s 1 -d T (C(s) × x s ) u 1 (11)
where C(s) is the path curvature vector. As a consequence, a Lyapunov-stable controller was proposed [START_REF] Dahroug | Visual servoing controller for time-invariant 3D path following with remote centre of motion constraint[END_REF] as:

u 1 = αx s + βd (12)
where α depicts the tangential velocity and β is a control gain for correcting the lateral error. Since u 1 is controlled by the surgeon, α and β are not independent from each other. It turned out that a good strategy is to make α dependent on β since the surgical constraint is to respect the path rather than the progression along it. As a consequence, α was chosen as:

α = u 1 2 -β 2 d 2 if β 2 d 2 < u 1 2 0 else (13)
With that, priority is given to the reduction of the lateral error and then, if possible, to progress along the path. Experiments on a real robot (Figure 9) achieved a path following error of some 300µm while maintaining the distance of the tool body to the incision point within 300µm, in a suboptimal hardware set-up. Through this paper, several examples of path following in the domain of microrobotics were reported. We tried to explain for all cases, the user requirements and why decoupling path following in a closed-loop control from the manual or openloop control of the velocity along the path was relevant. Short technical descriptions and some results were given, as an appetizer for reading the full papers and converting to path following in microrobotics.

  This work was supported by µRALP (FP7 ICT project no. 288663), Labex ACTION (ANR-11-LABX-01-01), MultiFlag (ANR-16-CE33-0019-02) and µRoCS (ANR-17-CE19-xxxx-xx).

Fig. 1 .

 1 Fig.1. Geometric description of the laser spot steering[START_REF] Séon | Decoupling path following and velocity profile in vision-guided laser steering[END_REF]

Fig. 2 .

 2 Fig.2. Endoscopic Laser Surgery set-up[START_REF] Séon | Decoupling path following and velocity profile in vision-guided laser steering[END_REF]

Fig. 9 .

 9 Fig. 9. Set-up for path following under single port access constraint


Benjamin Perret
email: benjamin.perret@esiee.fr

Jean Cousty

Silvio Jamil

Ferzoli Guimarães

Yukiko Kenmochi

Laurent Najman

Removing non-significant regions in hierarchical clustering and segmentation

First algorithm for removing non-relevant regions in hierarchies of partitions.

• Efficient algorithm with O(n log(n)) time complexity.

• The regions of the simplified hierarchy are regions or union of regions of the initial hierarchies.

• Evaluation on natural image analysis and illustrations on hierarchical clustering of data points.

• The method can be used as a pre-or post-processing step to enhance the quality of hierarchical segmentation algorithms.

Introduction

Many algorithms for image segmentation or data clustering contain a step that removes unimportant regions or clusters. In this paper, we are dealing with the more general problem of removing unimportant regions from a hierarchy of partitions, while still preserving the hierarchical partition structure. This is a common problem that appears in many different situations. For example, constrained connectivity [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF] solves the chaining problem well-known as one of the issues with minimum spanning tree based approaches, but it may creates a series of small undesirable regions in situation where there is a ramp discontinuity (see [START_REF] Soille | Constrained connectivity and transition regions[END_REF] for an analysis of this particular case).

One way to achieve such a hierarchical simplification would be to extract all the possible segmentations from the hierarchy, and for each one of them, remove the non-important regions by merging these regions with one of their neighbours. One of the issues is that those merging steps have to be performed in a consistent way, so that the set of simplified segmentations is still a hierarchy. Another important issue is that such a process would be slow.

In the literature on transformations of hierarchical segmentations [START_REF] Guigues | Scale-sets image analysis[END_REF][START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF][START_REF] Cardelino | A contrario selection of optimal partitions for image segmentation[END_REF][START_REF] Bosilj | Hierarchical image representation simplification driven by region complexity[END_REF][START_REF] Xu | Hierarchical segmentation using tree-based shape spaces[END_REF], there is not guarantee that unimportant regions are removed from the hierarchy. For example, small regions (with small area) can appear at very high level in the hierarchical tree, and the methods do not remove them. Thus, there is a need for post-processing the hierarchy. To the best of our knowledge, no algorithm has ever been presented for performing such a task. In order to provide such an efficient algorithm for removing unimportant regions from a hierarchy of partitions, we rely on the framework proposed in [START_REF] Cousty | Hierarchical segmentations with graphs: Quasi-flat zones, minimum spanning trees, and saliency maps[END_REF], where the equivalence between various hierarchical representations (dendrograms, saliency maps or minimum spanning trees) is demonstrated (see Section 2). As shown in Fig. 1, our algorithm makes use of these different representations to efficiently achieve its goal. This algorithm has been briefly introduced in the appendix of [START_REF] Guimarães | Hierarchizing graph-based image segmentation algorithms relying on region dissimilarity: the case of the felzenszwalb-huttenlocher method[END_REF], but a detailed analysis and clear explanations were missing; they are provided in Section 3. Furthermore, an empirical evaluation demonstrating its practical effectiveness is performed in Section 4.

Basic notions for graph-based hierarchy processing

Any hierarchy can be equivalently represented by sets as series of nested partitions or with a characteristic function defined on the edges of a graph and called a saliency map. The core of the hierarchy simplification method which we propose in this article and which is precisely described in Section 3 considers the saliency map representations of the hierarchies. In this section, we provide the formal definitions of the set representation of hierarchies and of the saliency maps. We also highlight how one can switch between the set and the functional representations of a hierarchy.

Hierarchies of partitions

In this article, the symbol V denotes a finite set which stands for the working space. In applications to image analysis, it can be for instance the set of all image pixels or superpixels. A partition of V is a set P of nonempty disjoint subsets of V whose union is V. Each element of P is called a region of P.

A hierarchy on V is a sequence H = (P 0 , . . . , P ) of partitions of V such that, for any i ∈ {1, . . . , }, any region of the partition P i-1 is included in a region of P i .

Figure 2 (a) and (c) illustrate an example of a hierarchy and of its dendrogram, respectively. Dendrograms are commonly used in applications. Intuitively, the dendrogram of a hierarchy embeds the inclusion relationship between the regions of this hierarchy. More precisely, it is a tree where the nodes correspond to the regions of the hierarchy and where each region R is linked to the largest (non-empty) regions of the hierarchy which are proper subsets of R, called the children of R.

Saliency map

Any hierarchy can be represented by an edge-weighted graph [START_REF] Cousty | Hierarchical segmentations with graphs: Quasi-flat zones, minimum spanning trees, and saliency maps[END_REF] spanning the elements of the space V. We provide in this section the definition of such a representation called a saliency map. A graph (spanning V) is a pair G = (V, E) such that E is a subset of the set of all unorderd pairs of distinct elements of V, i.e., the set E is a subset of {{x, y}

⊆ V | x y}. If G = (V, E) is a graph, each element of V is called a vertex of G, and each element of E is called an edge of G. A subgraph of a graph G = (V, E) is a graph (V , E
) such that V and E are subsets of V and of E, respectively. If X is a graph, its vertex and edge sets are denoted by V(X) and E(X), respectively.

A sequence (x 0 , . . . , x ) of vertices of a graph G is called a path from x 0 to x in G if any two successive vertices in the sequence form an edge of G, i.e., for any i in {1, . . . , }, the unordered pair {x i-1 , x i } is an edge of G. A graph is connected whenever there is a path from any of its vertices to every other one. Let G be a graph, by extension, we say that a subset R of V is connected

(for G) if the subgraph of G induced by R is connected, i.e., the subgraph (R, {{x, y} ∈ E(G) | x ∈ R, y ∈ R}) is connected. A connected component of G is a subset R of V
which is connected and maximal for this property: any proper superset of R is not connected.

In the sequel of this article, we assume that the space V is structured by a graph G = (V, E). For instance, in applications to image analysis, if the set V contains the set of all pixels or superpixels of an image, the edge set E can be obtained by any pixel or superpixel adjacency relation such as the one induced by the classical 4-, 6-or 8-adjacency relations. Furthermore, we will also assume that any hierarchy on V is connected for G meaning that any region of any considered hierarchy is connected for the graph G. These assumptions correspond to the situations which are the most often encountered in hierarchical image analysis. However, they can be dropped by considering that the graph G is the complete graph on V so that any subset of V is always connected. In such case the notion of a saliency map, whose definition is recalled hereafter, corresponds exactly to the notion of ultrametric distance which is well known in classification [START_REF] Leclerc | Description combinatoire des ultramétriques[END_REF].

A map w from E into the set R of real numbers is called a weight map on G. For any edge u of E, the value w(u) is called the weight of u, and the pair (G, w) is called an edgeweighted graph. Given a graph G = (V, E) and a hierarchy H on V, we show below how to define the saliency map Φ G (H) of H from E to R, which is an equivalent representation of the hierarchy H; knowing H one can infer Φ G (H) and, conversely, knowing Φ G (H) one can recover H.

Let us consider a hierarchy H = (P 0 , . . . , P ) on V. The saliency map of H is the map Φ G (H) from E to L = {0, . . . , } ⊂ R, such that the weight of any edge u = {x, y} of G for Φ G (H) is the largest value λ in L such that x and y belong to two distinct regions of P λ . Figure 2 (b) shows the saliency map Φ G (H) of the hierarchy H given in Figure 2 (a).

There is a bijection between the set of all hierarchies on V and the set containing every map which is the saliency map of a hierarchy (see Theorem 1 of [START_REF] Cousty | Hierarchical segmentations with graphs: Quasi-flat zones, minimum spanning trees, and saliency maps[END_REF]). In the next section, we present the quasi-flat-zone transform, denoted by QF Z which is the inverse of Φ G and allows to recover the hierarchy H knowing only its saliency map Φ G (H). These two transforms, namely Φ G and QF Z, make it possible to treat a hierachy either in a "set-oriented domain" (left part of Figure 1) or in the "saliency map domain" (right part of Figure 1).

An algorithm for computing the saliency map of any hierarchy H in linear time with respect to the size of the graph G, i.e., O(|V| + |E|), is described in [START_REF] Cousty | Hierarchical segmentations with graphs: Quasi-flat zones, minimum spanning trees, and saliency maps[END_REF]. This algorithm can be sketched as follows:

1. preprocess H for least common ancestors searches; 2. for each edge u = {x, y} of G taken in any order, 2..1 find the least common ancestor R of x and of y in the dendrogram of H; 2..2 set the weight of u to the level of R in the hierarchy.

Quasi-flat zone hierarchy

Quasi-flat zone transform [START_REF] Nagao | Region extraction and shape analysis in aerial photographs[END_REF][START_REF] Meyer | Morphological scale-space representation with levelings[END_REF][START_REF] Cousty | Hierarchical segmentations with graphs: Quasi-flat zones, minimum spanning trees, and saliency maps[END_REF] maps any edgeweighted graph into a hierarchy. In particular, if the departing map is the saliency map of a hierarchy, this transform allows to recover the initial hierarchy. As the hierarchy simplification method which we propose in Section 3 treats the hierarchies from their saliency maps, the quasi-flat zone transform allows us to recover the hierarchy associated with the saliency maps produced by our simplification method (see the overview diagram of Figure 1). Intuitively, this transform considers the series of the connected component partitions induced by the successive level sets of the edge-weight map.

Given an edge-weighted graph (G, w) and a value λ ∈ R, the λ-level set of E for w is defined by w λ (E) = {u ∈ E | w(u) < λ} and its associated subgraph (V, w λ (E)), denoted by w λ (G), is called the λ-level graph of G for w. The set of all connected components of w λ (G), denoted by C(w λ (G)), is a partition of V called the λ-level partition of G for w. Given an edge-weighted graph (G, w), the quasi-flat zone hierarchy QF Z(G, w) of (G, w) is then the finite sequence of all λ-level partitions of G for w, ordered by increasing values of λ, namely,

QF Z(G, w) = (C(w λ (G)) | λ ∈ R).

Hierarchy simplification with an attribute criterion

The proposed simplification method is based on a regional attribute, such as region area (size) and contrast, which measures the significance of any region. It aims at transforming an initial hierarchy into a new one such that:

• the new hierarchy does not contain any region with an attribute value below a given threshold;

• the regions of the new hierarchy are either regions of the initial hierarchies or regions obtained by merging adjacent regions of the initial hierarchy.

As mentioned in the introduction, in order to efficiently perform such simplification, the hierarchies are represented by weight maps. More precisely, we consider the (spatially and functionally) minimal representation of a hierarchy introduced in [START_REF] Cousty | Hierarchical segmentations with graphs: Quasi-flat zones, minimum spanning trees, and saliency maps[END_REF]: it consists of a minimal weighted subgraph (in terms of inclusion relation on graphs) whose quasi-flat zone hierarchy is precisely the hierarchy that we aim to represent. Such minimal representation of a hierarchy can be obtained by considering first the graph G weighted by the saliency map of the given hierarchy and then restricting it to one of its minimum spanning trees (Theorem 12 in [START_REF] Cousty | Hierarchical segmentations with graphs: Quasi-flat zones, minimum spanning trees, and saliency maps[END_REF]), leading to a weighted tree (T, w). The core of the method is then to produce a new weight map w for this graph T , standing for the saliency map representing the resulting simplified hierarchy. In order to produce such map, the edges of the tree are considered in any order. For each edge {x, y}, the largest region of the hierarchy which contains x but not y and the largest region of the hierarchy which contains y but not x are analyzed. If the attribute value of one of these two regions is below the given threshold, then the two regions must be merged. This is done by setting to 0 the weight of {x, y} for w . On the contrary, if the attributes of both regions are above the given threshold, then the two regions must be kept and we replicate the weight of {x, y} for w into w . In order to efficiently implement this method, a fundamental operation consists of finding the largest region of a hierarchy which contains one extremity of an edge but not the other. This can be done with the help of a data structure called a binary partition tree by altitude ordering. Hence, before giving a precise presentation of the simplification algorithm in Section 3.2, we first present, in Section 3.1, binary partition trees by altitude ordering together with a simple algorithm to compute them.

Binary partition tree by altitude ordering

Binary partition trees by altitude ordering (BPTAOs) are deeply related to Kruskal's minimum spanning tree algorithm [START_REF] Kruskal | On the shortest spanning subtree of a graph and the traveling salesman problem[END_REF].

Function B.AddNode(x) 1 B.Parent[x] = -1; 2 B.Size += 1; Function B.FindRoot(x) 1 while B.Parent[x] ≥ 0 do x B.Parent[x]; 2 return x Function B.CreateParent(x,y) 1 i B.Size; /* index for the new node */ 2 B.AddNode(i); 3 B.Parent[x] i; 4 B.Parent[y] i; 5 B.LeftChild[i] x; 6 B.RightChild[i] y;
More precisely, the BPTAO data structure can be seen as the (tree-based representation of a) hierarchy of partitions of V obtained during Kruskal's minimum-spanning-tree algorithm. A formal definition of this structure can be found in [START_REF] Cousty | Constructive links between some morphological hierarchies on edge-weighted graphs[END_REF] and algorithms to construct them are presented in [START_REF] Najman | Playing with kruskal: Algorithms for morphological trees in edge-weighted graphs[END_REF]. In this article, for the sake of completeness, we present a simple algorithm to construct it. However, the reader interested into a more efficient construction is refered to [START_REF] Najman | Playing with kruskal: Algorithms for morphological trees in edge-weighted graphs[END_REF].

This simple construction of a BPTAO from an edge weighted graph (G, w) is given in Algorithm 1. It corresponds to a particular implementation of Kruskal's algorithm. The auxiliary functions called in Algorithm 1, namely AddNode, FindRoot and CreateParent, are also described below Algorithm 1. In Algorithm 1, we initially consider a partition into singletons (Line 2) which is the first level of the BPTAO. Then, when an edge is selected by Kruskal's algorithm, we build the next level by merging the largest regions containing the vertices of the selected edge {x, y} (Lines 3-9). In terms of tree, the newly created region R is a new node of the BPTAO B, which becomes the parent of the two nodes associated with the merged regions (Line 7). There is a direct relation between the newly created region R and the edge {x, y} that is considered for the merging which creates the region R. In Algorithm 1, we observe that the edge {x, y} is stored at an index e (see Line 8) of the array L MS T and that the index of the region R in the tree data structure B is n + e = |V| + e (Line 7), allowing to keep track of the relation between the edge {x, y} and the region R for further processing. When the algorithm terminates, we then obtain:

• a minimum spanning tree of (G, w) whose edges are stored, following a non-decreasing order of weight (called an altitude ordering), in the array L MS T ;

• a tree B, called the BPTAO of (G, w) associated with L MS T , whose non-leaf nodes correspond to the edges the minimum spanning tree produced by Kruskal's algorithm (Line 8) and whose leaves correspond to the vertices of G (Line 2);

• an implicit mapping between the nodes of the BPTAO B and the vertices and edges of the minimum spanning tree. Any node of B stored at an index between 0 and n -1 is mapped to the vertex of the graph G at the same index, whereas any node of B with an index i between n and 2n-1 is mapped to the edge of the minimum spanning tree stored in

L MS T [i -n].
Figure 3 illustrates the relationship between a minimum spanning tree of G and its associated BPTAO B.

It should be also noticed that, if the quasi-flat zone hierarchy QF Z(G, w) is a binary hierarchy (i.e., each region is either a singleton or the result from the merging of exactly two regions), then it is equal to the BPTAO produced by Algorithm 1 [START_REF] Cousty | Constructive links between some morphological hierarchies on edge-weighted graphs[END_REF]. Otherwise, the hierarchy QF Z(G, w) can be straightforwardly recovered from B as shown in [START_REF] Cousty | Constructive links between some morphological hierarchies on edge-weighted graphs[END_REF][START_REF] Najman | Playing with kruskal: Algorithms for morphological trees in edge-weighted graphs[END_REF].

For a more efficient implementation of Algorithm 1, readers are referred to [START_REF] Najman | Playing with kruskal: Algorithms for morphological trees in edge-weighted graphs[END_REF]. Provided that the edges of the graph G are either already sorted or can be sorted in linear time, the efficient algorithm of [START_REF] Najman | Playing with kruskal: Algorithms for morphological trees in edge-weighted graphs[END_REF]

O(|E(G)| × α(|V(G)|))
, where α is the extremely slowly growing inverse of the single-valued Ackermann function.

Hierarchy simplification algorithm

The hierarchy simplification method is precisely described, with the help of Playing with Kruskal algorithm (namely Algorithm 1), in Algorithm 2. In the first line of the algorithm, an MST of the saliency map w of a given hierarchy H and its associated BPTAO B are obtained from Algorithm 1. After calculating the attribute for every region R in B (Line 2), we can efficiently carry out the two main steps of the method for each edge u ∈ E(G), thanks to the two structures L MS T and B: (1) get the attribute values of the associated connected components in Lines 5 and 6, and (2) set the new edge weight w (u) depending on the verification of the attribute criterion for the two regions merged by the edge u (Lines 8 and 9). It can be observed that Line 7 uses the mapping between the nodes of B and the edges of the considered MST, which was presented at the end of Section 3.1 and which is illustrated with the green arrows in Fig. 3.

It can be observed that, as presented in Figure 1, the hierarchy given to Algorithm 2 and the one resulting from its execution are in the form of saliency maps denoted by w and w respectively. The tree-based representation of the simplified hierarchy resulting from Algorithm 2 can be obtained by computing the quasi-flat zone hierarchy of w for the MST stored in L MS T . Such computation can be done, for instance, with the algorithm presented in [START_REF] Najman | Playing with kruskal: Algorithms for morphological trees in edge-weighted graphs[END_REF].

Illustrations and assessments

This section presents qualitative and quantitative assessments of the proposed method. Our tests focus on two different hierarchical segmentation methods: the quasi-flat zone hierarchy (QFZ, see Section 2.3) and the watershed hierarchy by area (WS-Area). Watershed hierarchies were first proposed in [START_REF] Beucher | Watershed, hierarchical segmentation and waterfall algorithm[END_REF][START_REF] Najman | Geodesic saliency of watershed contours and hierarchical segmentation[END_REF][START_REF] Meyer | The dynamics of minima and contours[END_REF] and have since been formalized in the context of minimum spanning forests [START_REF] Cousty | Incremental algorithm for hierarchical minimum spanning forests and saliency of watershed cuts[END_REF]. Intuitively, the WS-Area hierarchy of an edge-weighted graph is obtained by sequentially filtering the edge weights of the graph with area closings of increasing sizes and then computing the sequence of watershed segmentations of these filtered edge weights.

Then, we consider two regional attributes to simplify those hierarchies:

1. the area of a region, defined as the number of vertices in the region; and 2. the frontier strength of a region, defined as the mean weight of the edges linking the region with its sibling, i.e., the edges on the frontier between the two children of the parent region.

The area attribute of each region can be computed in linear time from the BPTAO by traversing the tree from the leaves to the root, the area of the leaves being 1 and the area of a non-leaf node being the sum of the area of its two children. The frontier strength can also be computed in linear time, by traversing the edges of the graph G, finding the lowest common ancestor of the two vertices of the edge in the BPTAO (this query can be done in constant time thanks to a linear time pre-processing of the tree [START_REF] Bender | The lca problem revisited[END_REF]) and accumulating the edge weights in this region.

The area attribute is used to identify non significant nodes in the QFZ hierarchy, while the frontier strength attribute is used in conjunction with the WS-Area hierarchy.

We first present illustrations of non-significant node removal on hierarchies built on point clouds and images. Then, we present extensive quantitative assessments of the benefits of our procedure for natural image analysis.

Illustrations

We first demonstrate the effectiveness of the proposed method on the hierarchical analysis of two simulated 2D point clouds (see Figure 4) 1 . Each point cloud is generated from three random distributions corresponding to three classes. We then consider the graph induced by the Delaunay triangulation 1 All the illustrations presented in this section can be reproduced using the Python Notebooks available at https://higra.readthedocs.io/en/ latest/notebooks.html.

of the points and we weight the edges by the Euclidean distance between the points. We observe that very small regions are branching at very high levels in the dendrogram of the QFZ hierarchy. Hence, the partition containing three regions in the hierarchy fails to correctly recovers the three clusters. By removing non-significant nodes from the QFZ hierarchy based on an area attribute (nodes containing less than 7 points are considered non-significant), we ensure that the hierarchy does not contain any small region anymore (neither at high nor at low levels). We observe that the partitions containing three regions in the simplified hierarchy correctly recovers the three clusters.

Another illustration of the effectiveness of the proposed method on hierarchical natural image analysis is demonstrated in Figure 5. In this Figure, the saliency map of a hierarchy (Section 2.2) is represented in the 2D Khalimsky grid [START_REF] Arbelaez | Contour detection and hierarchical image segmentation[END_REF][START_REF] Najman | Geodesic saliency of watershed contours and hierarchical segmentation[END_REF]: in this representation, the brightness of a contour is inversely proportional to the number of partitions of the hierarchy it belongs to, i.e., dark contours are the strongest ones. We can observe that in the QFZ hierarchy, most strong contours represent very small regions located on thick transitions between different regions of the images. When the saliency map is plotted in the 2D Khalimsky grid, this suppression of small regions looks like a sharpening, in other words, thick and blurred transitions become sharp. On the contrary, we can see that WS-Area hierarchy already produces thin contours. However, it also produces a lot of non-significant contours in large homogeneous regions. After a region removal procedure with a small contour strength from the WS-Area hierarchy, most spurious contours disappear.

Quantitative assessment

This section presents a quantitative assessment of the proposed method on natural image analysis. We first explain the assessment methodology, the evaluation measures, and the image datasets. Then, we give the results comparing the QFZ and WS-Area hierarchies to their simplified counterparts. Finally, we also compare our results with the one obtained by the transformation of a hierarchy into its optimal cut hierarchy for the piecewise constant Mumford-Shah energy [START_REF] Guigues | Scale-sets image analysis[END_REF].

Methodology. We mainly follow the supervised assessment framework proposed in [START_REF] Perret | Evaluation of hierarchical watersheds[END_REF]. We give an overview of the quality measures and readers can refer to the provided references to get detailed descriptions. The assessment framework relies on three types of measures to encompass various aspects of hierarchical representations:

1. precision-recall and F-measure on boundaries (FB) [START_REF] Arbelaez | Contour detection and hierarchical image segmentation[END_REF].

This measure evaluates the quality of the boundaries of each partition of a hierarchy with respect to a groundtruth segmentation. To evaluate a hierarchical method on a whole dataset, two aggregated measures are then defined: 1) the optimal image scale (OIS) measuring the best achievable score when taking the optimal partition in each hierarchy, and 2) the optimal data-set scale (ODS) measuring the best achievable score when taking partitions at the same level (the optimal scale) in every hierarchy; 2. fragmentation curves on the bidirectional-consistency error (BCE) [START_REF] Perret | Evaluation of hierarchical watersheds[END_REF]. The fragmentation of a partition is defined Input image QFZ hierarchy simplified QFZ hierarchy WS-Area hierarchy simplified WS-Area hierarchy Fig. 5. Removal of non-significant nodes of the QFZ hierarchy and of the WS-Area hierarchy on 4 images of the BSDS 500 dataset [START_REF] Arbelaez | Contour detection and hierarchical image segmentation[END_REF]. For each image, we show from left to right: the input image, the saliency map of the QFZ hiearchy, the saliency map of the simplified QFZ hierarchy, the saliency map of the WS-Area hierarchy, and the saliency map of the simplified WS-Area hierarchy.

as the number of regions in the partition divided by the number of regions in the ground-truth. The fragmentation curve on BCE evaluates the quality of the regions of partitions of the hierarchy as the fragmentation increases, also with respect to a ground-truth segmentation. We consider two categories of partitions that can be extracted from a hierarchy: the partitions of the hierarchy (horizontal cuts), and the optimal partitions constructable from regions taken from any partition of the hierarchy (the optimal non-horizontal cuts). Two aggregated measures are defined: the area under the curve for optimal cuts (FOC) and the area under the curve for horizontal cuts (FHC);

3. object detection measure [START_REF] Perret | Evaluation of hierarchical watersheds[END_REF][START_REF] Perret | Evaluation of morphological hierarchies for supervised segmentation[END_REF]. This last measure is based on supervised object detection with markers (one marker for the object and one for the background) and tries to describe an object as a set of regions taken from any partition of the hierarchy. It quantifies how well a specific object of a scene can be retrieved with different levels of in-formation given on its position. Markers are automatically generated from the ground-truth and corresponds to: erosions of the ground-truth object/background masks (Er), skeletons of the ground-truth object/background masks (Sk), and the frame of the image (Fr). Three combination of background-foreground markers are considered: Er-Er, Fr-Sk, and Sk-Sk. The quality of a detection is measure with its Jaccard index.

The precision-recall curves and fragmentation curves are evaluated on the test set of the Pascal Context dataset [START_REF] Mottaghi | The role of context for object detection and semantic segmentation in the wild[END_REF] which consists of a pixel-wise segmentation of the last 2 498 images of the Pascal VOC '10 [10] validation set. The object detection measure is evaluated on the MS-COCO [START_REF] Lin | Microsoft COCO: Common Objects in Context[END_REF] dataset. Each object of the dataset is processed independently leading to a total of 291 875 objects from the 40 504 images of the MS-COCO 2014 validation set. The edge weights of the graph of an image are then defined as the mean gradient value of its two extremities, the gradient being obtained with the structured edge detector [START_REF] Dollár | Fast edge detection using structured forests[END_REF]. In order to evaluate the benefits of the proposed method we propose two comparisons:

1. QFZ hierarchy versus a simplified QFZ hierarchy where small regions have been removed. The area threshold is expressed as a fraction of the total number of pixels in the image. 2. WS-Area hierarchy versus a simplified WS-Area hierarchy where regions with a common weak frontier have been merged. The strength threshold assumes that gradient values are normalized between 0 and 1.

Table 1 shows the results obtained with QFZ hierarchies. We can see that the removing of small regions provides significant improvements for the three measures. A threshold level of 0.4% or 0.8% (between 150 and 400 pixels on the tested images) offers a good compromise on the different measures.

Table 2 shows the results obtained with WS-Area hierarchies. In this case, the results are more contrasted. While a suppression of weak contours can provide significant improvement on precision-recall curves and fragmentation curves, the effect can be rapidly detrimental to the object detection measure. This issue can be due to the fact that the MS-COCO dataset contains a lot of poorly resolved objects with weak contours that can be deleted by the proposed method. However, we still see that a hierarchy simplification with moderate threshold values (between 0.05 and 0.1) improves all the considered quality measures.

Finally, the last lines of Tables 1 and2 show the results obtained by the transformation of a hierarchy into its optimal cut hierarchy [START_REF] Guigues | Scale-sets image analysis[END_REF]. We recall that this transformation modify the level of the nodes of a hierarchy such that each partition of the transformed hierarchy is optimal for the piecewise constant Mumford-Shah energy whose regularization parameter is equal to the level of the partition. We can see that this transformation provides very good results on the WS-Area hierarchy where it can identify incorrect contours, thanks to the rich information provided by the Mumford-Shah energy, and then push them down to the bottom of the hierarchy. It is however unable to deal with the small regions present close to the top of the QFZ hierarchy as pushing them down to the bottom would require to completely collapse the hierarchy. Nevertheless, we can see that the combination of the two transformation methods, the proposed simplification strategy followed by the transformation into optimal cut hierarchy, on QFZ (last line of Table 1) gives the best result. This further support the idea that the proposed method can be used as a pre-or post-processing step to enhance the quality of hierarchical segmentation algorithms.

Conclusion and perspectives

In this paper, relying on the framework developed in [START_REF] Cousty | Hierarchical segmentations with graphs: Quasi-flat zones, minimum spanning trees, and saliency maps[END_REF], we have provided a generic solution to the common problem of removing non-significant regions from a hierarchy of partitions. The experiments demonstrate that applying this algorithm does indeed improve the results in a number of situations. Future work will combine our approach with probability functions (e.g., attention saliency) or some other criterion relying on deep learning techniques to achieve state-of-the-art results.

Fig. 1 .

 1 Fig. 1. A flowchart of the proposed method for removing non-significant regions from a given hierarchy H and obtaining a new hierarchy H .

Fig. 2 .

 2 Fig. 2. (a) An example of hierarchy of partitions, H = (P 0 , P 1 , P 2 , P 3 ) where P 0 = {{a}, {b}, {c}, {d}, {e}, { f }}, P 1 = {{a, b}, {c, f }, {d}, {e}}, P 2 = {{a, b, d}, {c, e, f }}, P 3 = {{a, b, c, d, e, f }}; (b) the saliency map of H on a graph G, Φ G (H); (c) the dendrogram of the hierarchy H.

Algorithm 1 : 6 if r x r y then 7 B 8 L 9 e

 16789 Playing with Kruskal Data: An edge-weighted graph ((V, E), w) Result: An array L MS T to store the edges of an MST of G in non-decreasing order of weight with respect to w Result: Its associated BPTAO B 1 e 0 ; /* Initialize the index for L MS T */ 2 foreach x ∈ V do B.AddNode(x) ; /* Assuming that V = {0, . . . , n -1} */ 3 foreach {x, y} ∈ E in non-decreasing order of w do 4 r x B.FindRoot(x); 5 r y B.FindRoot(y); .CreateParent(r x , r y ); MS T [e] {x, y}; += 1;

7 u 8 if a 1 ≥

 781 L MS T [n -|V|]; m and a 2 ≥ m then w (u) w(u); 9 else w (u) 0;

5 Fig. 3 .

 53 Fig.3. Given a weighted graph, its minimum spanning tree T (whose edges are thick and gray) is represented by the binary partition tree by altitude ordering B (in blue). Each leaf node corresponds to a vertex of T while each non-leaf node n i of B corresponds to an edge of T ; the correspondences are depicted in green arrows.

Fig. 4 .

 4 Fig. 4. Removal of non-significant nodes on the QFZ hierarchical clustering of two point clouds (first and second lines). For each graph, we show from left to right: the graph with the three ground-truth clusterings of the graph vertices (red, green, and blue), the dendrogram of the QFZ hierarchy, the clustering into 3 classes for this dendrogram, the dendrogram of the simplification of the QFZ hierarchy, and the clustering into 3 classes for this simplified dendrogram. Note that the colors used to represent clusterings are arbitrary and do not represent an explicit correspondence between two different clusterings.

  has a quasi-linear time complexity, Algorithm 2: Hierarchy simplification by attribute Data: A graph G = (V, E) that is the working space Data: The saliency map w of a hierarchy H Data: An attribute threshold value m Result: A (saliency) saliency map w defined on the edges of an MST T of (G, w)

	1 Calculate the ordered edge array L MS T and its associated
	BPTAO B (Algorithm 1);
	2 Calculate the region attribute of each node n of B and
	store it in A[n];
	3 foreach non-leaf node n of B /* n iterates over
		the set {|V|, . . . , 2|V| -1}	*/
	4 do
	5	a 1 A[B.LeftChild[n]];
	6	a 2 A[B.RightChild[n]];

Table 1 .

 1 Comparison between the QFZ hierarchy, the simplified QFZ hierarchy at a given threshold, and the hierarchy of optimal Mumford-Shah (MS) cuts of QFZ.

		ODS	FB	OIS	BCE FOC FHC Mean Median OD	Mean score
	QFZ 0.479 0.477 0.358 0.358 0.500	0.550	0.454
	QFZ + simplification 0.4%	0.525 0.580 0.510 0.464 0.552	0.613	0.541
	QFZ + simplification 0.8%	0.537 0.589 0.533 0.475 0.533	0.579	0.541
	QFZ + simplification 1.6%	0.517 0.602 0.550 0.483 0.505	0.524	0.530
	QFZ + simplification 3.2%	0.515 0.543 0.541 0.479 0.467	0.448	0.499
	QFZ + MS cuts [11] 0.525 0.523 0.368 0.368 0.503	0.551	0.473
	QFZ + Simplification 0.8% + MS cuts 0.595 0.637 0.548 0.505 0.528	0.569	0.564

Table 2 .

 2 Comparison between the WS-Area hierarchy, the simplified WS-Area hierarchy at a given threshold , and the hierarchy of optimal Mumford-Shah (MS) cuts of WS-Area.

	ODS	FB	OIS	BCE FOC FHC Mean Median OD	Mean score
	WS-Area 0.512 0.591 0.588 0.440 0.518	0.552	0.534
	WS-Area + simplification 0.05 0.522 0.592 0.589 0.445 0.519	0.554	0.536
	WS-Area + simplification 0.08 0.527 0.596 0.591 0.452 0.519	0.554	0.540
	WS-Area + simplification 0.10 0.530 0.599 0.591 0.457 0.518	0.553	0.541
	WS-Area + simplification 0.15 0.541 0.604 0.593 0.470 0.511	0.539	0.543
	WS-Area + simplification 0.20 0.541 0.605 0.592 0.482 0.490	0.503	0.536
	WS-Area + MS cuts [11] 0.535 0.585 0.615 0.514 0.531	0.576	0.559

Results. Each image of the test datasets presented in the previous section was first transformed into a 4-adjacency graph.

Acknowledgments

The authors are grateful to CNPq (Universal 421521/2016-3 and PQ 307062/2016-3), FAPEMIG (PPM-00006-16) and PUC Minas for the financial support to this work. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brasil (CAPES) -Finance Code 001 and by CAPES/COFECUB 88887.191730/2018-00.


Mathilde Negri
email: mathilde.negri@lirmm.fr

Maximilien Servajean
email: maximilien.servajean@lirmm.fr

Benjamin Deneu
email: benjamin.deneu@inria.fr

Alexis Joly
email: alexis.joly@inria.fr

Location-Based Plant Species Prediction Using A CNN Model Trained On Several Kingdoms -Best Method Of GeoLifeCLEF 2019 Challenge

This technical report describes the model that achieved the best performance of the GeoLifeCLEF challenge, the objective of which was to evaluate methods for plant species prediction based on their geographical location. Our method is based on an adaptation of the Inception v3 architecture initially dedicated to the classification of RGB images. We modified the input layer of this architecture so as to process the spatialized environmental tensors as images with 77 distinct channels. Using this architecture, we did train several models that mainly differed in the used training data and in the predicted output classes. One of the main objective, in particular, was to compare the performance of a model trained with plant occurrences only to that obtained with a model trained on all available occurrences, including the species of other kingdoms. Our results show that the global model performs consistently better than the plant-specific model. This suggests that the convolutional neural network is able to capture some inter-dependencies among all species and that this information significantly improves the generalisation capacity of the model for any species.

Introduction

Predicting a list of the most likely species present at a given location can be very useful. First, it could improve species identification processes and tools by reducing the list of candidate species observable at a given location (whether automated, semi-automated or based on classical field guides or flora). More generally, it could facilitate biodiversity inventories through the development of location-based recommendation services (typically on mobile phones) as well as the involvement of non-expert nature observers. Finally, it could be used for educational purposes through biodiversity discovery applications providing innovative features such as contextualized educational pathways.

This challenge is related to the problem of Species Distribution Modelling (SDM) in ecology. SDM goal is to predict the spatial distribution of a species over a territory [START_REF] Kearney | Mechanistic niche modelling: combining physiological and spatial data to predict species' ranges[END_REF], in our case we used spatial positions and also environmental data. SDM have become increasingly important in the last few decades for the study of biodiversity, the ecology of conservation, landscape management, preservation of rare and/or endangered species, measurement of human impact or climate change on species, etc.

Concretely, the objective of SDM is to infer the spatial distribution of a given species, and they are often based on a set of geolocalized occurrences of that species (collected by naturalists, field ecologists, nature observers, citizen sciences project, etc.).

Recently, SDM based on deep neural networks have begun to appear [START_REF] Botella | A deep learning approach to species distribution modelling[END_REF]. These experiments have shown that they can have a good predictive power, potentially better than the models used conventionally in ecology such as MaxEnt [START_REF] Phillips | Maximum entropy modeling of species geographic distributions[END_REF]. Deep neural networks can learn complex non-linear transformations in a wide variety of fields. In addition, they provide an opportunity to learn an area of environmental representation common to various species, which stabilizes predictions from one species to another and improves them globally. Finally, spatial patterns of environmental variables often contain useful information for species distribution, but are generally not considered in conventional models. On the contrary, convolutional neural networks effectively use this information and improve their predictions.

In this paper, we present a study to evaluate a convolutional neural network to determine the ecological preferences of species through ranges of environmental image patches provided as input (temperature, soil type, etc.) as part of the GeoLifeCLEF challenge.

-Section 2 gives an overview of the various data we used to build our model.

-Section 3 provides the detailed description of our model.

-Section 4 presents the results of the experiments and their analysis.

Data

Occurrences

The data set protocol is explained in the challenge protocol [START_REF] Botella | Overview of geolifeclef 2019: plant species prediction using environment and animal occurrences[END_REF].

plantnet data set: composed of 2,367,145 plant species occurrences with uncertain identifications, because they come from the automatic identification of plant images, from the Pl@ntnet application. glc 18: Global Biodiversity Information Facility (GBIF) data set (same data set as last year challenge), it is composed of 291,392 occurrences of 3,336 plant species observed on French territory between 1835 and 2017.

trusted data set: a sample of plantnet data set without uncertainty.

no plant: 10,618,839 species occurrences from the GBIF database, of other kingdoms (such as mammals, birds, amphibians, insects, fungus etc.). We have removed occurrences that didn't match any environmental rasters, for example birds in the middle of the sea.

Then, we filter the database to get only plant species that were given by the challenge.

We used different methods to structure our train data set from each data set.

full noplant: we train our model with all the data (plantnet, glc 18 and no plant) we had including animal observations and species from other kingdoms.

full prediction: we only used the Pl@ntnet data set and glc 18 data set.

filter predictions: we try a random weighted selection scheme, the same as the one explained in the challenge notes for the test set. For each occurrence s i in our data set, we compute a weight w i which corresponds to the importance of this plant over an area of 2 kilometres. Then we did a random sampling on the z i number between 0 and the maximum weight, and we kept the occurrence only if z i < w i . We did this sampling over 2 data set: the plantnet data set with uncertain identification set at 70%, we had a train set around 32,000 occurrences after sampling, and another with the data set pl trusted, we had around 27,000 occurrences after sampling.

Environmental rasters

We have associated 33 rasters windows cropped into each corresponding global environmental rasters provided by the challenge as input features to each occurrence. [START_REF] Botella | Overview of geolifeclef 2019: plant species prediction using environment and animal occurrences[END_REF] An extraction protocol is given by the challenge. 5These environmental rasters were constructed from various open data sets including Chelsea Climate, ESDB soil pedological data, Corine Land Cover 2012 soil occupation data, CGIAR-CSI evapotranspiration data, USGS Elevation data (Data available from the U.S. Geological Survey.) and BD Carthage hydrologic data. All these data try to best represent the environment where the plant is observed. We print them on a given area in the Figure 1, we can see the diversity of information we have for one occurrence near Montpellier.

In the following, we generally denote to x X an occurrence, each x being associated to a spatial position p(x) in the spatial domain D, a species label y(x) and an environmental tensor g(x) of size 64x64x33. We denote as P the set of all spatial positions p covered by X. It is important to note that a given spatial position p 0 ∈ P usually corresponds to several occurrences x j ∈ X, p(x j ) = p 0 observed at that location. In the training set, up to several hundreds of occurrences can be located in the same place.

The environmental data provided [START_REF] Botella | Overview of geolifeclef 2019: plant species prediction using environment and animal occurrences[END_REF] is composed of tensors of 64x64X33 pixels. The corresponding 64x64 pixel matrices can be processed as classical image channels provided as CNN input. Most of them are continuous variables such as average temperature, altitude or distance to water as we can see them in the Figure 1. Thus, the corresponding 64x64 pixel matrices can be processed as classical image channels provided as CNN inputs. Some of the variables are ordinal type (such as ESDB v2), they can be considered as additional channels of the CNN in the sense that the order of pixel values is not significant. For categorical variables like the Corine Land Cover soil type variable provided within GeoLifeCLEF. This variable can take up to 48 different categorical values, but 3 of them are not used therefore we have kept 45 categorical variables. The order of these values has no meaning. Consequently, we preferred unstacking the corresponding channel into 45 different binary channels, then concatenate it with the other 32 tensors. Thus, instead of having a tensor of 64x64X33 pixels, we do have a tensor of 64x64x77 pixels.

Convolutional Neural Network

It has been previously shown in [START_REF] Botella | A deep learning approach to species distribution modelling[END_REF] that Convolutional Neural Networks (CNN) may reach better predictive performance than classical models used in ecology.

Our model attempts to predict the most likely species to be observed based on environmental features learned. Our structure is very different from last year challenge [START_REF] Deneu | Location-based species recommendation using co-occurrences and environment-GeoLifeCLEF 2018 challenge[END_REF]. We used an Inception V3 convolutional neural network [START_REF] Szegedy | Rethinking the inception architecture for computer vision[END_REF], it is made up of symmetric and asymmetric building blocks, including convolutions, average pooling, max pooling, concatenates, dropouts, and fully connected layers. In addition to this model, we added a first layer of 77 layers, to match the size of the environmental raster and a dropout to avoid over fitting.

For the last layer we used a softmax and a cross entropy loss. The softmax layer is computed as:

p j = e hj K k=1 e h k (1)
where y k are the scores inferred by the net for each class, and K the number of different classes that means the number of species we trained our model with, the softmax activation for a class y j depends on all the scores in y.

Then we used a cross entropy loss [START_REF] De Boer | A tutorial on the crossentropy method[END_REF]:

Loss = - K k=1 y k log(p k ) (2)
with y k a binary indicator if class label k is the correct classification for the observation and p k the probability of the observation is of class k.

Learning set up and parameters: all our experiments were conducted using PyTorch deep learning framework 4 and were run on a single computing node equipped with 4 Nvidia GTX 1080 ti GPU. We used the Stochastic Gradient Descent optimization algorithm with a learning rate of 0.1, a momentum of 0.9 and a batch size of 300. We trained our model over multiple epochs, then to export our result we took the epoch where the model gets the highest Top 1 score in our validation set. For our best runs (full noplant predictions and full predictions) we trained only over 20 epochs due to the number of data and for the two other runs we trained our model over 180 epochs.

Result

For GeoLifeCLEF 2019 they decided that the main evaluation criteria will be the accuracy based on the first 30 answers, also called Top 30, i.e. the function scoring 1 when the right specie is in the 30 first answers, and 0 if not. And as a For this challenge we submitted 4 different runs. We can see the global result in Figure 2 and Table 1.

-The first one: full noplant predictions is a CNN trained with every data set (plantnet, no plant, glc 18). We used pl@ntnet database without any filters (2,367,145 plant species occurrences), train and test occurrences from last year's challenge (file GLC 2018.csv) and occurrences from other kingdoms (such as mammals, birds, amphibians, insects, fungus etc.) from the GBIF database (file noPlant.csv). There are 34,375 different classes, so this CNN has a larger last layer output, but for the test we let it choose only the plant classes which are provided by the challenge. This method showed the best result with a Top 30 of 0.1769 and a MRR of 0.031. -The second submission, full predictions is a CNN trained with plantnet data set and the data set of last year challenge. This means that we only have plants to train our model, we have 3,859 classes. It is the 5 th submission of the challenge, after SaraSi models with a top 30 of 0.1364. -The two other submissions were not as good (plcomplete predictions and inception glc19 filter predictions), we sampled our train set as explained in the challenge protocol for the test set. Therefore, we have reduced our data, our model didn't capture as much information as the other models.

In this case, it seems that the more data available to the CNN, the better the results. It seems that our model is learning from the species from another kingdoms that live in the same area. Indeed, our model train with more than only plant occurrences showed a better result for predicting plants. Here, the CNN has more classes but can classify plants even better. We can deduct from this that our CNN network not only captures environmental information but also, the interaction between different species, we can imagine that it can also learn from the information of the species that live around the plants and perhaps from the ecological niche of an occurrence.

Conclusion

This paper describes our participation in GeoLifeCLEF challenge to evaluate location-based species prediction models. We compared Convolutional Neural Network trained with different data sets. The main conclusion of our study is that the convolutional neural network model is the most efficient model and it can learn information not only from environmental rasters but also from interaction between species of other kingdoms. Indeed, it achieved the best performance of the whole GeoLifeCLEF challenge when we trained our model with other kingdom observations.

In future work, we will attempt to better understand what information the CNN does capture from different data sets and how it could be improved.

Fig. 1 .

 1 Fig. 1. Environmental rasters

Fig. 2 .

 2 Fig. 2. Top30 accuracy score of every participant of GeoLifeCLEF 2019

Table 1 .

 1 Table of Result of GeoLifeCLEF 2019 challenge

	Rank runId top30 runName	participant
	1	27007 0.1769 full noplant predictions	LIRMM
	2	27086 0.1687 RUN2 GRINNELL FULL INITIAL	SaraSi
	3	27087 0.1653 RUN3 GRINNELL TRANS	SaraSi
	4	27088 0.1646 RUN1 GRINNELL TESTSPECIES INITIAL	SaraSi
	5	27006 0.1364 full predictions	LIRMM
	6	26997 0.1342 submit xgb spatial 4x4 all	SSN CSE
	7	26996 0.1288 submit xgb spatial allnoclc	SSN CSE
	8	27013 0.1273 submit xgb dep3 1	SSN CSE
	9	27069 0.1268 submission sel 4x4	SSN CSE
	10	27012 0.1263 submit xgb 4x4 all dep3	SSN CSE
	11	27070 0.1227 submission sel 1	SSN CSE
	12	27064 0.1198 submission1x1	SSN CSE
	13	27067 0.1135 submission4x4	SSN CSE
	14	27124 0.1135 Lot Of Lof 2	Lot of Lof
	15	27089 0.1110 RUN4 ELTON TRANS	SaraSi
	16	27082 0.1090 RUN0 ELTON FULL INITIAL TESTSPECIES SaraSi
	17	26988 0.1063 submit xgb spatial	SSN CSE
	18	27123 0.0984 Lot Of Lof 3	Lot of Lof
	19	27063 0.0864 Lot Of Lof 1	Lot of Lof
	20	26875 0.0844 submission	SSN CSE
	21	27102 0.0834 rfspatial	SSN CSE
	22	26821 0.0570 submit	SSN CSE
	23	27004 0.0470 plcomplete predictions	LIRMM
	24	27005 0.0465 inception glc19 filter predictions	LIRMM
	25	26968 0.0205 run 14	sergiu atodiresei
	26	26964 0.0191 run 10	sergiu atodiresei
	27	26961 0.0190 run 7	sergiu atodiresei
	28	26971 0.0184 run 17	sergiu atodiresei
	29	26967 0.0180 run 13	sergiu atodiresei
	30	26960 0.0168 run 6	sergiu atodiresei
	31	27062 0.0159 run 20	sergiu atodiresei
	32	26958 0.0146 run 3	sergiu atodiresei
	33	26970 0.0102 run 16	sergiu atodiresei
	34	26969 0.0099 run 15	sergiu atodiresei
	35	26972 0.0089 run 18	sergiu atodiresei
	36	26963 0.0079 run 9	sergiu atodiresei
	37	26965 0.0068 run 11	sergiu atodiresei
	38	26926 0.0067 run 4	sergiu atodiresei
	39	26973 0.0064 run 19	sergiu atodiresei
	40	26959 0.0063 run 5	sergiu atodiresei
	41	26962 0.0062 run 8	sergiu atodiresei
	42	26957 0.0061 run 2	sergiu atodiresei
	43	26966 0.0058 run 12	sergiu atodiresei
	44	26956 0.0058 run 1	sergiu atodiresei

http://www.github.com/maximiliense/GLC19


Pauline Guyot

El-Hadi Djermoune

Thierry Bastogne

Bruno Chenuel
email: b.chenuel@chru-nancy.fr

A Single ECG Lead-Based Oscillation Index for the Quantification of Periodic Breathing in Severe Heart Failure Patients

Keywords: Periodic breathing, Sinusoidal model, Cheyne-Stokes respiration, Electrocardiogram, Heart failure

Periodic breathing is a sleep-disordered breathing characterized by the alternation of central hypopneas/apneas and hyperventilation, and is associated with increased mortality in patients with severe heart failure in most studies. In this paper, we present a new strategy to detect mild to severe patterns of periodic breathing using a single electrocardiogram signal in patients with severe heart failure. We first compute three time series, extracted from the ECG signal namely Heart-Rate Variability, R-Wave Amplitude and Mean Cardiac Axis. Then, these series are used to estimate an oscillation index that can quantify periodic breathing through time and a one-minute decision is made using an experimental thresholding to decide whether periodic breathing is absent or present. Eight patients with normal to severe periodic breathing are selected to test our method. The results obtained are compared to those performed by experts.

I. INTRODUCTION

In patients with severe chronic heart failure, periodic breathing (PB) during sleep or daytime has been found to be a strong predictor of poor prognosis [START_REF] Lorenzi-Filho | Cheyne-Stokes respiration in patients with congestive heart failure: causes and consequences[END_REF]. Contrary to obstructive sleep apnea-hypopnea (OSA) characterized by a partial or complete obstruction of the upper airway during sleep, central sleep apnea-hypopnea (CSA) including PB, is an instability of the ventilation control leading to diminished or absent respiratory effort. It may be resulting from prolonged circulation delay or increased chemoreceptor sensitivity [START_REF] Francis | Quantitative general theory for periodic breathing in chronic heart failure and its clinical implications[END_REF]. PB can be separated into two classes, ventilation without apnea and ventilation with apnea called Cheyne-Stokes respiration (CSR) [START_REF] Lorenzi-Filho | Cheyne-Stokes respiration in patients with congestive heart failure: causes and consequences[END_REF], [START_REF] Francis | Quantitative general theory for periodic breathing in chronic heart failure and its clinical implications[END_REF]. Both classes are characterized by a crescendo-decrescendo pattern of ventilation, alternating hyperventilation and central hypopneas/apneas. Overnight multichannel polysomnography is commonly used for diagnosis and for the identification of patients with PB [START_REF] Rovere | Clinical relevance of short-term day-time breathing disorders in chronic heart failure patients[END_REF]. However, it is an expensive and timeconsuming test that is not available in all medical structures, and when it is, those units are typically overloaded. Simplified alternatives such as in-home respiratory polygraphy have been widely used and accepted even though its accuracy is uncertain [4]. Other methods have been proposed to reliably detect respiratory efforts, among them several are based on the electrocardiogram (ECG). One rests on the beat-to-beat variations of R wave-to-R wave (RR) intervals due to respiratory sinus arrhythmia, a phenomenon indexed on the cardiac vagal function that shortens RR intervals during inspiration and increases it during expiration [START_REF] Yasuma | Respiratory sinus arrhythmia: why does the heartbeat synchronize with respiratory rhythm?[END_REF]. A second approach called ECG-derived respiration relies on the change of distance between ECG electrodes and heart during respiration. It affects the amplitude of QRS complex and the mean cardiac electrical axis [START_REF] Moody | Derivation of respiratory signals from multi-lead ECGs[END_REF].

Typically, diagnosis for PB is established by computing the apnea-hypopnea index (AHI, number of apneas or hypopneas per hour). However, it only provides quantitative information on the respiratory signal and lacks any qualitative aspect such as the nature of events or the amplitude of the modulation when central. Some methods have been proposed to quantify the amplitude of the oscillations on a respiratory signal based on a standard amplitude modulation (AM) scheme with filters [START_REF] Tellez | eAMI: A Qualitative Quantification of Periodic Breathing Based on Amplitude of Oscillations[END_REF] or on a spectral decomposition algorithm of the instantaneous minute ventilation [START_REF] Pinna | Assessing the severity and improving the understanding of sleep-related breathing disorders in heart failure patients[END_REF]. Those methods are quite sensitive to noise and quantify the amplitude of the modulation but does not specify the characteristic of the oscillation like its instantaneous frequency, thus cannot confirm very typical pattern such as CSR. Although, all these methods are based on respiratory or thoraco-abdominal motion signals, none propose the quantification using electrocardiogram signals.

The purpose of this paper is to develop a new method able to detect, from a single ECG lead, mild to severe patterns of periodic breathing. Our main contribution relies on a sequential estimation of an oscillation index used as a biomarker to quantify the periodic breathing severity extracted from an ECG signal. This is a non-invasive and easy-to-use test. The remainder of the paper is organized as follows. The estimation problem is formulated in Section II. Then, the proposed algorithm for the estimation of repiratory series and the caculation of the oscillation index are detailed in Section III. In Section IV, the experimental results obtained on a panel of eight patients with chronic heart failure are compared with the opinions of sleep experts to assess the performances of the classification based on the oscillation index. Finally, the conclusions are drawn in Section V. There is a close relationship between the respiratory and the cardiovasculary systems, so when an instability such as a central apnea is introduced, the whole loop is disturbed. Instead of using a nasal cannula that can be noisy and unpleasant for the patient, our goal is to restitute and quantify the oscillation present in the respiratory efforts from a single electrocardiogram signal. Even if this test is an easy and noninvasive one, it is still an indirect respiration-recording and requires several processing steps to detect periodic breathing. Considering a single ECG lead, the respiration signal is computed then an oscillation index and its instantaneous frequency are estimated informing on the potential oscillation present in the signal due to PB. This oscillation index is considered as an efficient biomarker to graduate mild to severe patterns.

The proposed computational method can be decomposed into four successive steps: (i) detection of R-peaks from a single-ECG lead; (ii) computing of Heart-Rate Variability (HRV) and ECG-derived respiration series (RWA: R-wave amplitude, and MCA: mean cardiac axis); (iii) estimation of the oscillation index; and (iv) classification of each one-minute segment into normal or PB. These steps are sketched on Fig. 1 and are thoroughly described in the next section.

III. PROPOSED ALGORITHM

A. Estimation of respiratory series

The extraction of respiratory series from an ECG signal involves the detection of QRS complexes. This task is the easiest one in the processing of ECG signal and has been widely discussed in the litterature with methods based on quadratic energy [START_REF] Pan | A real-time QRS detection algorithm[END_REF], [START_REF] Hamilton | Quantitative investigation of QRS detection rules using the MIT/BIH arrhythmia database[END_REF], wavelet transform [START_REF] Merah | R-peaks detection based on stationary wavelet transform[END_REF], or adaptive filter [START_REF] Jain | QRS detection using adaptive filters: A comparative study[END_REF] to mention a few. Let x ECG (t) be the original ECG record. The algorithm we used for R-peaks detection is based on a sparse decomposition using a suitable redundant dictionary and can be summarized as follows.

1) Using a sliding window of 10 seconds, perform three steps: a) median filtering of 0.2 s for baseline correction [START_REF] De Chazal | Automatic classification of heartbeats using ECG morphology and heartbeat interval features[END_REF], b) high-pass filtering with a cut-off frequency of 50 Hz to remove low-frequency noise and other components of the ECG signal such as P and T waves, c) R-peak localization using Orthogonal Matching Pursuit (OMP) algorithm [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] with a dictionary composed of Coiflet-5 wavelet at level 3. The ventilation signal is shown for comparison purposes; it is not used in the processing method.

2) Elimination of multiple detection: when a peak is detected during the refractory period of 200 ms, it is disregarded. 3) Search back for missed R-peaks: if no R-peak is detected within 150% of the last RR interval, a search back is applied using OMP in the signal centered around the potential missed R-peak. Once all R peaks (l R ) are detected, three time series are generated: (1) The HRV series s HRV corresponding to interval time between two successive R-peaks; (2) The RWA series s RWA by computing R-peak amplitude which requires the suppression of the baseline wander through median filtering;

(3) The MCA series s MCA by computing the mean cardiac axis. The latter can be extracted from a baseline-free ECG lead by computing the QRS complex area on a fixed-length window of 80 ms around the detected R-peak. All signals are evenly resampled with spline interpolation. An example of these series is given on Fig. 2. In the sequel, the generic notation s x (t) will represent one of three series s HRV (t), s RWA (t) or s MCA (t).

B. Computation of the oscillation index

The next step is to characterize a potential oscillation present in the three extracted series to detect typical periodic breathing events. Unlike respiratory signals (e.g. see Fig. 2(a)), those series does not present a typical AM signal but the addition of two components as PB superimposes its lowfrequency oscillation on the cardiovascular system [START_REF] Leung | Influence of Cheyne-Stokes Respiration on Cardiovascular Oscillations in Heart Failure[END_REF]: the respiration component which is assumed to be a sinusoidal signal x r (t) = A r cos(2πf r t) whose frequency f r goes from 0.25 Hz to 0.33 Hz (from 15 to 20 respirations per minute), and the oscillation component which is either constant for a normal respiration or oscillating for a periodic breathing. Similarly, the latter is also assumed to be a sinusoidal signal x m (t) = A m cos(2πf m t + φ m ) whose frequency f m goes from 8 mHz to 30 mHz (a typical cycle of PB lasts from 30 s to 2 min). Note that f r f m . Finally, our model for the estimated series can be expressed as:

s x (t) = x r (t) + kx m (t) + α = A r cos(2πf r t) + kA m cos(2πf m t + φ m ) + α (1)
where α, k ∈ R are two constants. Let us denote by h m = kA m the oscillation index, a key parameter which graduates the amplitude level of the periodic breathing. The objective now is to estimate the oscillation index from s x (t) together with f m and φ m . To discard the respiration component, the signal is low-pass filtered using a cut-off frequency of 0.1 Hz. The output of the filter can be written in the following form:

s x (t) = kx m (t) + α (2) = a 0 e j2πf0t + a m e j2πfmt + a * m e -j2πfmt (3)
with f 0 = 0, a 0 = α, and a m = hm 2 e j2φ m , where j = √ -1 and a * m denotes the complex conjugate of a m . The phase φ m accounts for the delay introduced by the lowpass filter; this parameter is not of special interest here. Hence, our model is composed of three complex exponentials. Its parameters can be estimated by any subspace-based method. In this work, we use the Matrix Pencil method [START_REF] Hua | Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise[END_REF], [START_REF] Sarkar | Using the matrix pencil method to estimate the parameters of a sum of complex exponentials[END_REF] over a sliding window to observe the parameter variations. The window size t w has to be small enough for the stationarity assumption (the sinusoidal model) to hold. It is set to 2 minutes and the overlapping ratio ρ w between two successive windows is 80%. The oscillation index ĥm = 2|â m | and frequency fm are then computed and interpolated to get two continuous signals that can quantify and qualify through time the oscillation of the respiratory signal. This will be illustrated in section IV-B. t w is a critical parameter, to properly detect all events it must be set at least at 2 minutes and difference was found for t w ∈ [START_REF] Francis | Quantitative general theory for periodic breathing in chronic heart failure and its clinical implications[END_REF]4]. ρ w is also important to have a fine detection but a compromise must be established between the precision and the computation time.

C. Time-localization of periodic breathing

The next step is to decide whether the respiration is oscillating or not. Using a ROC analysis (Fig. 3) on the first two classes of the classification given by the two sleep experts (i.e. periodic breathing and non-periodic breathing), a threshold h 0 was used to detect a modulation of breathing sufficiently present to be pathological in the three series. The values of h 0 are 0.02, 0.03 and 0.032 for HRV, RWA and MCA, respectively. In parallel, fm has to be contained in an interval of frequencies [START_REF] Pinna | Assessing the severity and improving the understanding of sleep-related breathing disorders in heart failure patients[END_REF]30] mHz in which the oscillation component is associated to periodic breathing. If both ĥm and fm are classified as pathological for 1 minute, then a zone of PB (z osc ) is detected and the value of ĥm specifies the severity of the pathology.

IV. EXPERIMENTAL RESULTS

A. Data acquisition

This study is a retrospective analysis of data, which included adult patients referred to a sleep laboratory (University Hospital CHRU Nancy) for evaluation of suspected sleep disordered breathing. The study was approved by the Local Ethics Committee of the University Hospital of Nancy and informed consent was obtained from all subjects before they commenced participation. Subjects were seated comfortably on a chair in a quiet room, in a condition of relaxed wakefulness for about 30 minutes of recording. Inspiratory and expiratory flows were measured, the respiratory gas was continuously sampled from the pneumotachograph for the measurement of expired CO 2 and O 2 partial pressure. Electrocardiogram, oxygen saturation, thoracic belt respiration and blood pressure were also simultaneously recorded. The study involved a group of eight patients all presenting severe heart failure (with a left ventrical fraction ejection under 30%). Two tasks were asked to the sleep experts: the first one was to classify each minute of the ventilation signal based only on the visual aspect into three categories: (1) Non-PB, (2) PB and (3) Erratic breathing possibly PB. The second task was to establish diagnosis following international guidelines using all available signals. Two patients had severe PB with apnea (CSR), one patient exhibited a periodic breathing and five patients were classified as normal breathing. No patient was on opioids.

B. Results and discussion

Combining the value of the oscillation index ĥm and frequency fm allows to precisely describe the ECG-extracted respiratory signal. If ĥm is above the pathological threshold h 0 and if fm is contained in the interval mentionned before, then a pathological oscillation is sufficiently present to be marked as periodic breathing. However, it is important to combine both parameters together, if ĥm is above h 0 but fm is out of the interval or is not stable within, there is no oscillation zone marked. On the contrary, if fm is contained within the interval but ĥm less than h 0 , there is no readable oscillation in the envelope. The right way to analyse the parameters is:

1) to read carefully the values of ĥm . If it is steadily above h 0 then an oscillation is present, but we cannot conclude about its nature. If it is under h 0 , there is no oscillation in the signal; 2) to read fm signal: if it is continuously contained in the pathological interval, then an oscillation zone is detected; if fm is unstable (both within and outside the interval no more than 1 minute) or it is clearly outside the interval, then no oscillation is detected. In Fig. 2(a), the ventilation of a patient with severe PB is showed as a reference signal to compare with the three computed time series. Visual inspection allows us to differentiate oscillations associated with periodic breathing and breathing component in the computed series. It is also important to highlight that the pathological oscillation present in the ventilation signal, i.e. its envelope, can also be found in the Fig. 3. ROC analysis of the three time series to identify the optimal threshold to detect periodic breathing. HRV: Heart-Rate Variability series, RWA: R-peak Wave Amplitude series and MCA: Mean Cardiac Axis series series. Fig. 4(b)-(d) presents the HRV series and the estimated values ĥm and fm for a patient with Cheyne-Stokes respiration lowered or suppressed by the inhalation of CO 2 -enriched air (3% at t = 790 s and 6% at t = 1900 s) [START_REF] Lorenzi-Filho | Effects of inhaled carbon dioxide and oxygen on Cheyne-Stokes respiration in patients with heart failure[END_REF]. This record shows different levels of PB and the good performance of the algorithm with the different oscillation zones z osc detected by our method. ĥm is above the oscillation threshold and fm is included in the pathological interval when the patient presents CSR and ĥm and/or fm get out of the pathological intervals when the patient inhales a CO 2 -enriched air, known to be a respiratory regulator. The two parameters enable to conclude that both indices permit to correctly and locally detect and quantify the periodic breathing for HRV series. Compared to other methods that quantifies periodic breathing in ventilation signal [START_REF] Pinna | Assessing the severity and improving the understanding of sleep-related breathing disorders in heart failure patients[END_REF], our detection brings more information on the pattern of the oscillation as it can specify its instantaneous frequency, very typical and recognizable in the case of PB and CSR.

Using optimal thresholds presented in section III-C, Table I provides the numerical results obtained using the three series on all patients. MCA series achieve better overall performance than HRV and RWA for the oscillation detection on 1-minute window. When using all available signals, on the 31 minutes classified as erratic breathing possibly PB by the two experts, 21 minutes were clarified as PB and 10 as non-PB, mostly erratic breathing non resembling to PB. 23 minutes were classified as PB by the method: 21 minutes were correctly classified as periodic breathing in the MCA series and 2 minutes were false positives. Our method can help experts to detect PB events when erratic breathing hides a mild PB, difficult to detect with visual analysis.

In order to increase patients' comfort and to optimize medical care, methods to reduce the number of sensors have been developed and the use of ECG to monitor respiration have proven to be efficient [START_REF] Mazzanti | Validation of an ECG-derived respiration monitoring method[END_REF], [START_REF] Moody | Clinical validation of ecg derived respiration technique[END_REF]. In the case of heart failure patient, our method proposes to fetch different information from a single ECG signal, cardiac activity and sleep-disorder breathing, association strongly linked to mortality [START_REF] Khayat | Sleep disordered breathing and post-discharge mortality in patients with acute heart failure[END_REF]. The proposed detection method of PB has shown reliable results through our oscillation index estimation applied to patients with severe heart failure. Modeling the ventilation extracted from ECG with oscillation estimation leads to characterize its morphological aspect and provides an efficient numerical biomarker of PB severity. The higher h m , the more severe is the pathology. This parameter allows to locate patient's symptoms on the pathology continuum, to precisely locate abnormal events, and to graduate each epoch according to the value of h m instead of using the AHI, which is a single indicator for the whole signal.

V. CONCLUSION

We presented a new computational and sequential method to estimate an oscillation index, used as an efficient biomarker to detect and estimate severity levels of periodic breathing from a single ECG signal. All the components of the proposed method have been tested on a panel of 8 patients. The classification results showed promising performances of the proposed solution and demonstrated the proof of concept. A short-term perspective will focus on the application of the proposed solution on a larger panel of patients with other sleep-breathing disorders such as obstructive events. With a larger panel, it will also be possible to tune the threshold used on the oscillation index through, for example, cross-validation.

Fig. 2 .

 2 Fig. 2. Portion of record for a patient presenting Cheyne-Stokes Respiration. (a) Ventilation signal (b) Heart-Rate Variability series (c) R-peak Wave Amplitude series (d) Mean Cardiac Axis series (e) Electrocardiogram signal.The series HRV, RWA and MCA are computed from the ECG signal in (e). The ventilation signal is shown for comparison purposes; it is not used in the processing method.

TABLE I COMPARISON

 I OF NUMERICAL RESULTS FOR THE DETECTION OF PERIODIC

			BREATHING
		HRV	RWA	MCA
	Specificity (%)	84.00	75.63	83.27
	Sensitivity (%)	77.17	82.60	81.52
	Accuracy (%)	82.28	77.38	82.83

ACKNOWLEDGMENT

The authors wish to thank the virtual FEDER hospital project of Lorraine, France (FEDER-Etat-Région Hôpital Virtuel de Lorraine) for their valuable support.